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1.  INTRODUCTION 

The  work  done  under  the  contract  F44620-71-C-0049 
"Research  in  Seismology"  during  the  period  1 April  1971 
31  May  1975  can  be  grouped  into  three  broad  categories. 

1.  Seismic  source  mechanisms. 

2.  Earth  structure. 

3.  Seismic  wave  propgation,  scattering  and  attenuation. 

Under  each  one  of  the  three  areas  both  theoretical 

studies  and  applications  to  field  data  were  carried  out. 

In  source  mechanism  studies,  radiation  of  seismic  waves 
from  an  expanding  fault  model  was  calculated  numerically, 
and  source  mechanisms  and  focal  depths  of  earthquakes  in 
central  Asia  were  determined  from  surface  wave  spectra.  In 
earth  structure,  techniques  were  developed  to  invert  the 
travel  time  in  terms  of  the  three  dimensional  earth  models. 
These  techniques  were  utilized  to  determine  regional  struc- 
ture (under  LASA  and  under  California)  and  lateral  velocity 

variations  in  the  earth's  mantle. 

In  summarizing  research  efforts,  generally  the  results 
are  described  very  briefly.  When  the  material  is  available 
in  open  literature  in  the  form  of  published  papers,  refer- 
able theses,  or  papers  in  press  nearing  publication,  only 
the  abstracts  are  given.  Recently  completed  papers  are 
included  in  full  in  the  report.  In  addition  to  these,  a 
list  of  publications  and  a list  of  theses  completed  during 
the  contract  period  are  given  at  the  end. 
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SEISMIC  SOURCE  MECHANISMS:  THEORY 


2.1  nvnamics  of  an  expanding  circular  fault  by  R.  Madariaga 

abstract 

We  study  a plane  circular  model  of  a frictional  fault  using 
numerical  methods.  The  model  is  dynamic  since  we  specify  the 
effective  stress  at  the  fault.  In  one  model  we  assume  that  the 
fault  appears  instantaneously  in  the  medium;  in  another,  that  the 
rupture  nucleates  at  the  center  and  that  rupture  proceeds  at 
constant  subsonic  velocity  until  it  suddenly  stops.  The  total 
source  slip  is  larger  at  the  center  and  the  rise  time  is  also 
longer  at  the  center  of  the  fault.  The  dynamic  slip  overshoots 
the  static  slip  by  15%-35%.  As  a consequence,  the  stress  drop 
is  larger  than  the  effective  stress  and  the  apparent  stress  is 

less  than  one  half  the  effective  stress. 

The  far-field  radiation  is  discussed  in  detail.  We  distinguish 

three  spectral  regions.  First,  the  usual  constant  low  frequency 
level.  Second,  an  intermediate  region  controlled  by  the  fault 
size  and,  finally,  the  high  frequency  asymptote.  The  central 
region  includes  the  corner  frequency  and  is  quite  complicated. 

The  corner  frequency  is  shown  to  be  inversely  proportional  to  the 
width  of  the  far-field  displacement  pulse  which,  in  turn,  is 
related  to  the  time  lag  between  the  stopping  phases.  The  average 
corner  frequency  of  S waves  v®  is  related  to  the  final  source 
radius,  a,  by  v*  = .21  3/a.  The  corner  frequency  of  P waves  is 
larger  than  v?  by  an  average  factor  of  1.5. 
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2.2  Timp-Dependent  Strain  Accumulation  and__Release at  Island 

Arcs_i  implications  for  the  1946  Nankaido  EarthguaKe  by 

A.T.  Smith,  Jr. 

ABSTRACT 

Underthrusting  and  the  elastic-rebound  theory  are 
consistent  with  the  gross  static  deformations  after  earth- 
quakes for  island-arc  regions  such  as  Japan,  Alaska,  and 
Chile.  Yet  anomalous,  time-dependent  post-earthquake 
adjustments  suggest  additional  processes.  Here  the  astheno- 
sphere  or  mantle  becomes  the  element  that  both  determines 
the  post-seismic  deformation  and  controls  the  accumulation 
of  strain.  The  lithosphere  and  asthenosphere  represent  a 
coupled  system.  A large  earthquake  strains  the  entire 
system;  stress  relaxation  in  the  viscous  asthenosphere 
follows  and  allows  the  post-seismic  readjustments. 

The  convergence  zone  is  first  considered  as  a semi 
infinite,  elastic  plate  overlying  a viscoelastic  foundation. 
Analytic  solutions  for  short-period  deformations  place 
bounds  on  the  behavior  for  the  surface  deformations,  boundary 
conditions  on  the  fault  interface,  and  stress-propagation 
following  an  earthquake.  Detailed  models  are  then  considered 
using  a novel,  time-dependent,  finite  element  solution  for 
the  convergence  zor.e.  The  solution  avoids  propagation  o 
errors  in  time  and  readily  accommodates  inversion  theory. 

The  method  clearly  defines  the  behavior  for  realistic  models. 


Thus,  scaling  with  the  fault  depth  and  lithospheric  thick- 
ness controls  the  shape  for  simple,  planar  fault  models,  while 
the  time  scale  depends  on  the  asthenospheric  viscosity. 
Different  boundary  conditions  imposed  on  the  fault,  whether 
a constant  dislocation  or  a constant  stress  with  time, 
strongly  affect  the  resulting  deformations  and  stresses. 
Finally,  stresses  introduced  by  thermal  density  anomalies 
within  tha  descending  lithosphere  are  compared  to 
models  and  focal  mechanisms  near  Hokkaido. 

Generalized-matrix  inversion  theory  now  places  bounds 

on  the  effective  viscosity  and  fault  parameters  using 

geodetic  data,  focal  mechanisms,  and  tectonic  setting  for 

the  1946  Nankaido  earthquake  (M  8.2)  in  southwest  Japan. 

Using  the  assumption  of  stress  relaxation  in  the  astheno- 

sphere,  the  data  constrains  the  fault  geometry  to  a shallow 

15°  dip,  followed  by  a 60°  dip  from  26  km  depth  to  the  base 

of  the  lithosphere  at  60  km  depth.  Near  the  hi  oocenter 

the  slip  is  3 meters,  while  the  maximum  slip  is  less  than 

15  meters.  Other  models  with  constant  dip  or  shallower  dip 

beyond  25  km  do  not  satisfy  these  constraints.  The 
2 20 
viscosity  o{  the  asthenosphere  now  becomes  10  poise.  The 

results  suggest  segmentation  of  the  lithosphere  and 

deformations  that  generate  the  embayed  shoreline,  sedimentary 

basins  off  southwest  Japan,  seismicity,  and  focal  mechanisms. 
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2 . 3 A_Disl ocation  Approach  to  Plate  Interaction  by 

R.L.  Brown,  Jr. 

ABSTRACT 

A dislocation  can  be  described  in  terms  of  a surface  of 
discontinuity  or  the  line  which  circumscribes  this  surface. 

„e  have  applied  the  solutions  of  Yoffe  (1960)  and  Comninou 
(1973)  for  an  angular  dislocation  line  to  the  problem  of 
calculating  the  fields  due  to  general  polygonal  dislocations. 

Next,  a numerical  method  has  been  developed  explicitly 
for  finite  sources  (Finite  Source  Method  or  FSM)  which  allows 
the  computation  of  fields  from  a dislocation  that  penetrates 
several  layers  of  a layered  half-space.  The  speed  of  the  FSM 
allows  the  calculation  of  many  models  which  are  not  econom- 
ically possible  by  other  means.  It  is  used  here  to  model 
earthquakes  in  layered  media  and  plate  bottom  effects  due 

to  the  interaction  of  lithospheric  plates. 

Finally,  the  problem  of  the  mutual  interaction  of 
lithospheric  plates  in  relative  motion  has  been  posed  in 
terms  of  dislocation  theory  (anti-dislocations) . Dislocation 
models  of  various  portions  of  the  San  Andreas  fault  in 
California  are  proposed  and  evaluated  by  comparing  them  with 
seismic  and  geodetic  data.  We  find,  for  example,  that  fault 
creep  near  Hollister  acts  to  obscure  any  locking  at  depth 
and  that  as  much  as  70%  of  the  fault  could  be  locked  (down 
to  20  km)  and  still  be  consistent  with  the  geodetic  data. 


II 


The  models  also  suggest  that  the  depth  of  locking  (or 
non-slipping  portion  of  the  fault)  varies  from  10  to  80  kn 


along  the  San  Ai.’reas.  Under  San  Francisco  the  depth  of 


. _ on  4-n  4f)  km  while  iust  north  and  south 
locking  appears  to  be  20  to  su  Km  wnnc 


of  this  region  the  locking  is  from  10  to  15  km  deep.  Our 


models  are  also  indicative  of  a more  northerly  component  of 


motion  for  the  Pacific  with  respect  to  the  American  plate 


than  would  be  expected  if  the  San  Andreas  were  a simple 


strike-slip  fault.  South  of  Cholame  the  depth  of  locking 


begins  a rapid  increase  and  appears  to  lock  to  80  km  in  the 


Tejon  bend  portion  of  the  San  Andreas.  We  are  not  able, 
however,  to  distinguish  between  an  actual  locking  of  the 


fault,  capable  of  taking  high  stresses,  or  simply  a low 


stress  state. 


2 . 4 on  the?  Interaction  of  Two  Scales  of  Convection_  in_the 


Mantle  by  F.M.  Richter  and  B.  Parsons 


ABSTRACT 


A system  in  which  convection  takes  place  in  the  upper 
mantle  on  two  distinct  horizontal  length  scales  is  proposed, 


This  is  consistent  with  the  existence  of  the  plates  them- 


selves, the  relatively  constant  heat  flux  background  in 


older  ocean  basins,  and  the  knowledge  of  convection  in 


fluid  layers  gained  from  laboratory  and  numerical  experiments, 


h 
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The  large-scale  circulation  consists  of  the  plates  themselves 
and  the  return  flow  necessary  to  conserve  mass.  The  small- 
scale  flow,  analogous  to  Rayleigh-Benard  convection  or 
variants  of  this,  which  have  been  the  main  target  of  numerical 
study,  provides  the  necessary  vertical  heat  transport  in  the 
upper  mantle  that  supplies  the  required  heat  flux  at  the  base 
of  the  lithosphere.  The  depth  of  convection  is  taken  to  be 
down  to  the  650-km  seismic  discontinuity,  and  this  depth 
characterizes  the  horizontal  length  scale  of  the  small-scale 
convection.  This  system  is  studied  by  means  of  a set  of 
laboratory  experiments  that  explore  the  interaction  of  the 
small-scale  convection  with  the  large-scale  flow.  The 
experiments  show  the  plausibility  of  convection  on  two  scales. 
Furthermore,  they  suggest  that  beneath  fast-moving  plates 
(absolute  velocities  around  10  cm  yr  ) the  small-scale 
convection  will  align  itself  as  rolls  in  the  direction  of 
the  large-scale  flow  in  geologically  short  times.  However, 
beneath  very  slow  moving  plates  the  times  required  for  the 
alignment  of  convective  rolls  are  long  in  comparison  with 
times  over  which  no  changes  in  plate  motions  are  to  be  ex- 
pected. Here  the  convective  planform  is  more  likely  to  take 
the  form  of  upwelling  and  downwelling  spouts.  Thus  this 
simple  system  of  a convecting  layer  beneath  a moving 
boundary  contains  the  possibility  of  explaining  a wide  variety 

Observational  tests  of  the  consequences 


of  surface  features. 


* 


A 


of  the  two-scale  idea  are  suggested,  and  the  assumptions  on 
which  this  idea  is  based  are  critically  discussed. 


3.  SEISMIC  SOURCE  MECHANISMS:  ASIA 

3.1  static  and  Dynamic  Fault  Parameters  of  the  Saitama 

Earthquake  of  July  1/  1968  bV  K*  Abe 

The  source  mechanism  of  the  Saitama  earthquake  (36.07‘N, 
1 3 9 . 4 0 ° E , Ms  = 5.4)  of  July  1,  1968,  is  studied  on  the  basis 
of  P-wave  first  motion,  aftershock,  long-period  surface-wave 
data  and  low-magnification  long-period  seismograms  recorded 
in  the  near-field.  A precise  location  of  the  aftershocks 
is  made  using  P and  S-P  time  data  obtained  by  a micro- 
earthquake observatory  network.  The  synthetic  near-field 
seismograms  based  on  the  Haskell  model  are  directly  compared 
with  the  observed  near-field  seismograms  for  wave  and  ampli- 
tude to  determine  the  dynamic  fault  parameters.  The  results 
obtained  are  as  follows:  source  geometry,  reverse  dip  slip 

with  considerable  right-lateral  strike-slip  component:  dip^ 
direction,  KS'E:  dip  angle  30":  fault  dimension,  10  x 6 km  ; 
rupture  velocity,  3.4  km/sec  in  the  direction  S30“E:  average 
dislocation,  92  cm:  average  dislocation  velocity,  92cm/sec: 
seismic  moment,  1.9  • 1025  dyn-cm:  stress  drop,  100  bar. 

The  effective  stress  is  about  the  same  as  the  stress  drop. 


\ 
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For  major  earthquakes  in  the  Japanese  Islands,  the  disloca- 
tion velocity,  <D> , is  found  to  be  proportional  to  the  stress 
drop,  a.  This  relaticn  can  be  expressed  by  <D>  = 
where  B is  the  shear  velocity  and  y is  the  rigidity.  This 
result  has  an  importai.ce  in  engineering  seismology  because 
the  stress  drop  scales  the  seismic  motion  in  the  vicinity 
of  an  earthquake  fault. 


3 . 2 A Fault  Model  for  the  Niigata  Earthquake  of  1964  by 


K.  Abe 
ABSTRACT 

The  source  mechanism  of  the  Niigata  earthquake  of  1964 

(M  ^ 7.5)  is  studied  in  detail  on  the  basis  of  the  P wave 
s 

first  motion,  S wave  polarization  angle,  long-period  surface 
wave,  aftershock,  precise  levelling,  tide  gage,  and  tilt 
measurement  data.  The  long-period  surface  wave  data  and 
geodetic  data  are  interpreted  consistently  in  terms  of  a 
thrust  fault  reaching  the  earth's  surface  and  having  a dip 
of  56°  toward  N81°W,  a dimension  of  80  km  (length)  x 30  km 
(width),  and  an  average  dislocation  of  3.3  m.  The  seismic 
moment  is  3.2  x 1027  dyne-cm.  The  stress  drop  is  estimated 
to  be  70  bars.  This  value  is  not  very  different  from  the 
stress  drops  obtained  for  moderate  to  large  shallow  earth- 
quakes which  occurred  in  the  Japan  Islands.  The  fault  plane 


t 


geometry  obtained  here  is  slightly  different  from  that  deter- 
mined from  the  P wave  first  motions.  Combining  this  result 
with  the  weak  beginning  of  the  initial  P waves,  we  may  inter- 
pret the  entire  faulting  process  in  terms  of  a multiple 
faulting  which  consists  of  the  two  events:  the  initial 

localized  rupture  is  followed,  after  about  4 sec,  by  the 
major  faulting,  which  is  responsible  for  the  excitation  of 
long-period  waves. 

3 . 3 Thermal  and  Mechanical  Models  of  Conl:inent-Continent_ 

Convergence  Zones  by  P.  Bird  and  M.N.  Toksflz 
ABSTRACT 

The  thermal  regimes  of  continent-continent  convergence 
zones  are  modelled  by  a finite-difference  technique,  assuming 
that  there  is  some  subduction  of  continental  crust.  Gravity 
and  heat  flow  profiles  are  generated  from  the  thermal  models. 
Subducted  crust  and  slab  remain  cool  except  at  the  upper 
surface  where  frictional  heating  is  important.  Crustal  rocks 
may  be  metamorphosed  or  melted  by  friction  while  the  plates 
are  converging,  or  by  radioactive  self-heating  more  than 
30  m.y . after  the  plates  have  stopped.  In  the  former  case, 
melting  requires  frictional  shear  stresses  between  500  and 
2000  bars  at  a plate  velocity  of  5 cm/year.  At  lower 
velocities  the  upper  limit  of  frictional  stress  is  greater. 


The  total  work  performed  in  continental  subduction  may  be 
minimized  if  the  angle  of  underthrusting  becomes  more  shallow, 

changing  the  location  of  subduction. 

A model  for  the  geometry  of  oceanic  and  continental 
slabs  in  the  Zagros  Mountains  is  presented  which  satisfies 
gravity,  heat  flow,  seismic,  and  geological  constraints. 
Continental  underthrusting  is  beginning  along  a new  fracture 
at  the  edge  of  the  Persian  Gulf,  isolatina  the  Arabian 
continental  shelf  from  the  subduction  process.  Slippage 
along  this  fault  is  Pleistocene  and  probably  does  not  exceed 
30  km.  Subduction  of  continental  crust  at  the  Crush  Zone 
is  probably  insignif ioant . 


I 4 
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Introduction 

Subduction  of  continental  crust  as  the  cause  of  orogenies 
and  mountain-building  is  a concept  with  a long  history.  The 
geologic  evidence  is  clearest  in  the  Alpine-Himalayan  belt, 
and  as  early  as  1916  Emile  Argand  developed  a theory  of 
Alpine  orogeny  in  which  Africa  moved  North  and  overrode  Europe. 
Unlike  the  previous  geosynclinal  theory  of  Hall  (1860)  and 
Dana  (1873),  this  provided  an  explanation  for  the  asymmetry  of 
many  mountain  ranges,  in  which  deformation  proceeds  from  the 
mafic  to  the  silicic  side  with  a series  of  overthrusts  in  one 
direction.  More  recently,  the  theory  of  plate  tectonics  has 
suggested  a mechanism  and  allowed  us  to  identify  certain 
eugeosynclines  as  converging  plate  margins  towards  which  an 
opposing  continental  shelf  (the  miogeosyncline)  is  transported 
by  sea-floor  spreading  and  subduction  (Dewey  and  J.  Bird,  1970). 
One  natural  explanation  of  the  thick  low-density  roots  of 
mountain  ranges  is  that  the  collision  of  continents  is  followed 
by  some  amount  of  underthrusting  of  continental  crust. 

In  continental  convergence  belts,  particularly  the 
Himalayas  and  the  Zagros,  geology  records  both  the  oceanic 
and  continental  phases  of  subduction.  On  the  Northern  margin 
of  both  ranges  there  is  an  elevated  plateau,  intruded  by 
Voluminous  andesites  and  basalts  associated  with  the  former 
subduction  of  an  oceanic  plate  beneath  the  plateau.  Adjacent 
to  the  plateau  there  is  a suture  zone  with  steep  contacts, 
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deformed  turbidities  with  clastic  blocks  (trench  sediments) 
and  ultramafic  ophiolite  complexes  indicating  the  location  of 
the  first  contact  between  continents.  These  suture  zones, 
however, are  not  very  active  seismically  or  geombrphically . 
Instead,  the  elevation  of  the  mountain  range  occurs  in  a folded 
and  overthrust  zone  frcm  100  to  300  km  wide  on  the  other  side  of 
the  suture  zone.  These  deformed  zones  are  not  regionally  meta- 
morphosed, but  the  Himalayan  orogenic  belt  contains  scattered 
granites  which  may  originate  along  the  thrust  plane  from 
frictional  melting  as  the  crust  of  the  mountain  range  thickens 
by  under thrusting.  The  convergence  belts  are  terminated  by  a 
downwarping  alluvial  basin,  formed  over  the  advancing  plate  as 
it  bends  down  to  begin  subducting  under  the  range.  A possible 
sequence  of  development  of  these  features  is  given  by  Figure  1, 
in  which  the  Zagros  stage  is  shown  as  B and  the  Himalayan  stage 


as  C. 


The  location  of  crustal  subduction  is  further  marked  by 
seismicity  and  gravity  anomalies.  Earthquakes  occur  in  diffuse 
bands  dipping  away  from  the  margin  of  the  alluvial  foredeeps. 
Fault-plane  solutions  for  these  events  (Fitch,  1970;  Nowroozi, 
1972)  confirm  the  mechanisms  to  be  shallow  underthrusting.  The 
Bouguer  gravity  anomaly  fields  (Wilcox  et  al ♦ , 1972;  Qureshy, 
et  al.,  1974)  reveal  large  minima  elongated  Along  the  ranges 
and  centered  between  foredeep  and  suture  zone.  This  is 
additional  evidence  of  continental  subduction. 


1 
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In  this  paper  we  have  assumed  the  validity  of  plate 
tectonics  and  convergence  by  subduction  as  a starting  point. 

This  allows  calculation  of  the  critical  geophysical  parameters 
of  temperature,  density,  seismic  velocity,  and  stress,  and  it 
is  only  with  knowledge  of  these  fields  that  we  can  assess  the 
self-consistency  of  the  plate-tectonic  interpretation  of 
mountain-building.  With  this  goal,  we  have  considered  a number 
of  idealized  cases  illustrating  general  principles,  and  progressed 
to  an  interpretation  of  the  Zagros  Mountains,  in  which  the 
structures  described  are  just  beginning  to  develop. 

We  begin  with  a calculation  of  temperature  fields  in 
continental  subduction.  In  the  second  section  we  investigate 
the  effects  of  the  low  density  of  continental  crust  and  great 
thickness  of  shield  lithosphere,  and  the  ways  in  which  these  may 
force  the  location  and/or  the  angle  of  subduction  to  change. 

Finally,  we  present  a model  for  the  Zagros  Mts.,  where  the 
transition  from  oceanic  to  continental  subduction  has  been 
occurring  over  the  last  million  years. 

Calculation  of  Thermal  Models 

Thermal  modeling  of  downgoing  slabs  of  oceanic  lithosphere 
has  already  been  performed  by  Minear  and  Toksflz  (1970a)  and 
ToksGz  et  al.  (1971,  1973).  Using  the  same  technique,  we  have 
concentrated  on  the  details  of  the  shallow  structure  of  continental 


ft* 


-15- 


subduction  zones  where  the  vertical  displacement  is  assumed 
to  be  less  than  100  km.  We  begin  with  simple  hypothetical 
models  to  show  the  range  of  possible  thermal  regimes. 

To  calculate  the  thermal  models  of  converging  plates  of 
continental  lithosphere  we  solve  the  heat  equation 


C |?  = V (KVT)  + H 
9 1 


where  T = temperature,  t = time,  C is  heat  capacity,  K is 
thermal  conductivity,  and  H is  heat  production.  In  the  finite 
difference  solution  of  (1)  we  use  a 41  by  41  grid  with  a vertical 
increment  of  5 km.  The  details  of  the  computational  technique 
and  its  accuracy  have  been  extensively  discussed  by  Mine.ar  and 
Toksttz  (1970a,  1970b,  1971)  and  Toksttz  et  al.  (1971,  1973). 

The  parameters  of  primary  interest  in  these  calculation^ 
are  subduction  rate,  subduction  angle,  and  the  amount  of  shear- 
strain  heating.  Less  critical  parameters  are  held  constant 
from  model  to  model;  all  the  constants  used  are  presented  in 

Table  1 . 

Low  angles  of  subduction,  not  more  than  30°  in  dip,  were 
employed  because  fault-plane  solutions  in  continental  convergence 
sites  such  as  the  Himalayas  and  Zagros  have  consistently  shown 
shallow  underthrusting  (Fitch,  1970;  Nowroozi,  1972).  Also, 

4-he  great  buoyancy  of  continental  crust  will  tend  to  prevent 
subduction  at  large  dips,  as  will  be  discussed  below. 


Frictional  heating  (or  shear-strain  heating)  at  the  top 
of  the  slab  is  an  important  parameter.  These  two  terms  are 
used  interchangeably  in  this  paper  because  the  actual  mechanism 
is  unimportant  in  thermal  modelling.  Whether  there  is  stick- 
slip  behavior  or  viscous  creep  (Turcotte  and  Oxburgh,  1969)  the 
heat  production  per  unit  area  of  fault  is  just 

Q = fe t dZ  - tV  (2) 

where  Q is  in  ergs/cm^-sec,  t is  shear  stress  in  dynes/cm^, 
e strain  rate  and  V is  the  relative  plate  velocity  in  cm/sec. 

Although  t is  an  uncertain  parameter,  its  upper  and  lower 
bounds  can  be  established.  Stress  drops  in  earthquakes  provide 
a lower  limit.  Evidence  of  metamorphism  or  melting  in  each 
particular  range  can  be  used  to  fix  the  shear  stresses,  and  it 
will  be  shown  that  there  is  a universal  upper  limit.  A set  of 
theoretical  thermal  regimes  are  calculated  for  continental 
convergence  models.  The  results  strongly  depend  on  model 
parameters,  which  cannot  be  fixed  uniquely.  In  the  models 
discussed  below,  the  relative  effects  of  different  parameters 
are  illustrated.  A summary  of  models  and  model  parameters  are 
listed  in  Table  1. 

The  first  case  considered  is  that  of  slow  convergence  at 
1 cm/yr  and  at  the  relatively  large  dip  angle  of  30°.  Because 


continental  convergence  occurs  at  the  site  of  former  oceans 
(as  shown  in  Fig.  1)  rather  than  along  new  lineations  within 
a plate,  there  is  a previously-subducted  oceanic  slab  ahead 
of  the  continental  material.  This  oceanic  slab  absorbs  heat 
through  its  cold  upper  surface  while  subducting  and  creates  a 
cool  zone  in  the  overriding  lithosphere  before  the  continents 
come  together. 

The  thermal  regime  after  80  km  of  convergence  is  shown 
in  Figure  2.  Even  at  this  low  velocity  the  radioactivity  of 
the  crustal  material  has  little  effect.  The  frictional  heating 
rate  is  proportional  to  slab  velocity  (equation  2)  and  hence 
is  also  low.  Thus  the  continental  crust  remains  cold  except 
near  its  leading  surface  which  is  locally  heated  to  400 °C  and 
could  undergo  blueschist  metamorphism.  The  dominant  influx  of 
cold  material  is  also  reflected  in  the  heat  flow  above  the  slab, 
which  is  depressed  by  some  0.6  HFU  at  40  km  behind  the 
subduction  zone. 

A Bouguer  gravity  anomaly  has  also  been  calculated  for 
this  model  by  assuming  that  the  structure  is  flat  and  two- 
dimensional.  The  temperature  anomalies  are  converted  to  density 
anomalies  using  thermal  expansion  coefficients  of  3.0  x 10  5/°C 
for  the  upper  crust,  2.5  x 10  ^ for  the  lower  crust,  and 
3.2  x 10"5  for  the  mantle  (Skinner,  1966).  Because  the  oceanic 
slab  consists  mostly  of  cold  mantle  material,  this  is  the  only 
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contribution  to  its  density  anomaly.  The  continental  crust 
introduces  a large  compositional  density  anomaly  by  displacing 
heavier  mantle  rock,  and  this  overshadows  the  positive  thermal 
density  anomaly.  By  balancing  the  negative  crustal  and  positive 
thernal  mass  anomalies,  we  calculate  that  this  amount  of  sub- 
ducted crust  could  be  balanced  by  a minimum  oceanic  slab  only 
25  km  longer  than  that  of  Fig.  2.  This  is  assuming  a mantle 
thermal  expansion  of  6 x 10“5/°C  (Sleep,  1975) ; for  a lower 
coefficient  the  slab  must  be  longer. 

A second  calculation  at  a lower  dip  is  shown  in  Figure  3. 
The  convergence  rate  of  6 cm/yr  is  closer  to  the  rates  computed 
for  the  Zagros  and  Himalayan  ranges  from  sea-floor  spreading 
poles  and  rates  (LePichon,  1968) . Stress  has  been  slightly 
increased  on  the  shallow  parts  of  the  fault  where  most  earth- 
quakes seem  to  occur;  it  is  1 ki lobar  for  depths  between  80  and 
40  km,  1.2  kb  for  35  to  25  km,  and  1.5  kb  from  20  to  10  km. 

Frictional  heating  is  6 to  9 times  greater  than  in  the 
previous  case,  and  produces  a marked  high-temperature  zone  along 
the  top  of  the  slab  around  the  fault  surface.  Temperatures  in 
this  hot  zone  are  above  the  solidus  of  water-saturated  granite 
and  are  about  50°C  below  the  dry-muscovite-granite  solidus  of 
Brown  and  Fyfe  (1970) . These  solidus  temperatures  would  be 
reached  if  subduction  were  to  continue. 

Heat  flow  and  gravity  anomalies  associated  with  this  model 
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are  also  shown  in  Fig.  3.  There  is  a small  high  heat  flow 


anomaly  and  broad  gravity  anomaly.  From  these,  the  position 
of  the  slab  is  only  generally  indicated.  If  local  isostasy 


were  operative  during  the  subduction  process,  the  region  some 


100  km  behind  the  trench  line  would  be  uplifted  to  3 km  average 


elevation. 


This  model  illustrates  the  possibility  that  frictional 


heating  could  produce  strong  metamorphism  or  partial  melting 


in  a thin  zone  near  the  fault.  But  it  does  not  explain  the 


regional  metamorphism  which  seems  to  occur  late  in  the  history 


of  each  orogenic  belt.  To  test  the  possibility  that  this 
orogeny  is  produced  by  radioactive  heating  of  thickened  crust, 


we  allowed  the  calculation  in  this  second  case  to  continue  for 


100  m.y . with  no  motion  of  the  slab.  After  the  influx  of  cold 


surface  rock  stops,  the  extra  radioactive  heat  production  in 
the  region  of  overlapped  crust  eventually  dissipates  the 
regional  cold  anomaly  and  creates  a heated  region. 


The  difficulty  with  this  type  of  calculation  is  that 


vertical  motions  will  not  cease  with  the  end  of  convergence. 

A doubled  layer  of  crust  will  produce  dramatic  isostatic  uplifts 


and  rapid  erosion,  leading  to  the  eventual  destruction  of  the 


crustal  bulge.  Thus  the  time  range  over  which  this  modelling 


remains  valid  depends  directly  on  the  rate  of  erosion. 


The  temperature  field  at  the  end  of  30  m.y.  with  no  erosion 


corrections  is  shown  in  Figure  4A.  Dashed  lines  show  the  former 


outlines  of  the  slab.  The  solidus  line  at  the  base  of  the  plate 

' 

has  moved  upward  while  the  hot  region  in  the  base  of  the 

l ] 

overriding  lithosphere  has  d \scipated.  Thus  the  oceanic  slab 
begins  to  lose  its  identity.  An  elevation  of  temperatures  is 

r 

seen  in  the  slab  at  the  end  of  the  crustal  portion,  but  on  tho 
whole  the  region  is  still  anomalously  cold.  The  l'twer  portion 

of  the  crust  is  heated  above  6tO°C  and  severely  metamorphosed, 

•I 

but  nowhere  in  the  crust  is  the  solidus  exceeded.  Thus  the 

extensive  granite  bodies  and  high  grade  metamorphic  belts  such 

. 

as  those  contained  in  the  Urals  and  the  Piedmont  province  and 
Merrimack  Synclinorium  of  the  Appalachians  are  not  explained  by 
this  model.  They  could  be  produced  if  the  erosion  rate  were 
low  (or  the  radioactivity  higher). 

In  Figure  4B  the  same  model  is  shown  after  a quiescent 
period  of  100  m.y.  Here  the  inherited  cold  anomaly  has  disappeared 
and  there  is  a general  upwarping  of  isotherms.  Where  the  000°C 
isotherm  moves  into  the  base  of  the  subducted  crust,  general 
melting  and  formation  of  plutons  could  occur.  The  small  positive 
heat  flow  anomaly  of  0.2  KFU  which  is  predicted  would  be  augmented 
by  the  slow  upward  mass  transport  in  response  to  erosion  and  the 
rapid  mass  transport  of  rising  volatiles  from  the  region  of  melt. 

The  essential  problem  with  this  mechanism  for  metamorphism 
and  orogeny  is  the  excessive  time  required  to  raise  temperatures 

{ 

by  radioactivity  alone.  Naylor  (1971)  has  placed  a 30  m.y. 
maximum  duration  on  the  Acadian  orogeny  (in  eastern  Vermont) 


including  the  emplacement  of  granites,  and  if  the  orogeny 
resulted  from  continental  convergence  and  the  granites  from 
radioactive  heating  then  this  is  faster  than  predicted. 
Similarly,  Oxburgh  and  Turcotte  (1974)  found  that  simple 
overthrusting  and  radioactive  heating  were  insufficient  to 
produce  the  East  Alpine  metamorphisr  Thus  we  must  reserve 
the  hypothesis  that  tha  heat  of  orogeny  is  transported  to  the 
crust  from  below  rather  than  created  in  it. 

If  the  oceanic  slab  of  Figure  2 were  to  become  detached 
near  the  leading  edge  of  the  subducted  crust  (the  130  km  point) , 
it  would  continue  to  sink.  The  asthenosphere  at  over  1100°C 
would  rise,  at  a rate  comparable  to  plate  velocities,  to  take 
its  place.  Losing  heat  mainly  by  adiabatic  decompression,  it 
could  rise  to  only  50  km  depth,  diffuse  heat  to  the  crust  in 
under  lo  m.y.,  and  cause  a major  thermal  event. 


Mechanics  of  Continental  Convergence 

Before  it  is  possible  to  ascertain  the  maximum  possible 
extent  of  continental  subduction,  it  will  be  necessary  to  more 
accurately  quantify  the  various  forces  involved.  Body  forces 
may  be  quite  accurately  estimated  by  a combination  of  thermal 
modelling,  measurement  of  thermal  expansions,  and  observation 
of  surface  vertical  tectonics  (Smith  and  Toksflz,  1972).  Elastic 
and  viscous  bending  stresses  in  the  lithosphere  have  already 
been  estimated  by  observations  on  geologically  transient  loading 


t 


by  icecaps  (Walcott,  1970a),  glacial  lakes  (Walcott,  1970b), 
and  oceanic  subduction  (Smith,  1974;  Watts  and  Taiwan!,  1974). 

One  of  the  most  uncertain  factors  remaining  is  the  level  of 
frictional  stress  on  the  long  fault  surface  where  the  greatest 
strains  are  concentrated.  It  is  difficult  to  apply  laboratory 
experiments  to  this  problem  because  of  the  enormous  complexities 
of  continental  surface  rocks  and  the  uncertainty  as  to  how  much 
soft  sediment  and  water  are  entrained  into  the  fault.  One 
strength  of  thermal  modelling  is  that  it  can  set  definite  upper 
limits  on  the  possible  stresses  at  each  depth  in  this  complex  zone. 

The  thickness  of  the  fault  zone  is  not  critical  to  the 
thermal  history  as  long  as  it  is  much  less  than  the  scale  distance 
for  heat  conduction  during  the  subduction  phase.  Then,  the 
greatest  part  of  the  frictional  heat  produced  will  be  transferred 
to  the  slab  and  overlying  plate  rather  than  expended  in  heating 
the  shear  zone  itself.  Thus  the  temperatures  on  the  fault  will 
depend  only  on  the  gross  parameters  of  stress  and  plate  velocity, 
dip  angle,  and  initial  geotherm,  and  not  in  any  way  upon  the 
complex  mechanics  of  the  shear  strain  itself.  Since  all  these 
factors  except  stress  can  be  determined-  within  fractional  errors, 
what  remains  is  a strong  relation  between  shear  stress  and 
temperature. 

The  maximum  possible  shear  stress  at  any  depth  along  the 
fault  will  be  that  which  causes  the  temperature  to  reach  the 
solidus.  As  soon  as  partial  melting  of  the  rock  sets  in, 
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there  is  seif-regulating  effect  which  will  maintain  that 
temperature  by  reducing  the  effective  viscosity  or  coefficient 
of  frictin.  The  interesting  point  in  the  case  of  continental 
subduction  is  that  half  of  the  rocks  in  the  fault  zone  will 
have  a grossly  granitic  composition  and  hence  a very  low- 
temperature  solidus.  The  solidus  of  water-saturated  granite, 
which  is  flat  with  pressure  beyond  7 kb  (Tuttle  and  Bowen, 

1958)  is  not  particularly  relevant  because  significant  aim.unts 
of  water  (2%  by  weight)  would  be  consumed  by  the  formation  of 
partial  melts.  However,  the  solidus  of  Brown  and  Fyfe  (1970) 
for  dry  granite  containing  muscovite  may  be  more  appropriate 
for  the  convergence  zones.  As  shown  in  Figure  5,  this  solidus 
would  never  be  reached  along  the  fault  at  depths  of  less  than 
100  km  without  frictional  heating.  This  is  because  the  influx 
of  cold  material  depresses  the  fault-zone  geotherm  well  below 

the  base  geotherm. 

To  determine  the  shear  stresses  necessary  for  melting  we 
used  a combination  of  finite-difference  and  analytic  techniques. 
First,  we  used  the  slab  temperature  program  to  solve  for  temperatures 
on  the  fault  in  the  artificial  case  of  no  shear-strain  heating. 

This  includes  the  effects  of  slab  geometry,  non-linear  geotherm, 
adiabatic  compression,  etc.  Then,  since  the  finite-difference 
grid  is  too  course  to  accurately  represent  a thin  heat  source, 
we  used  asymptotic  analytic  formulas  to  calculate  the  values  of 
q (Z ) that  would  raise  these  fault  temperatures  to  the  solidus. 
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Figure  5 shows  the  temperatures  for  the  artificial  zero- 


stress  case  after  400  km  of  subduction.  This  distance  is  the 


approximate  width  of  the  continental-rise  sediments  on  an 


Atlantic-type  continental  margin  (Emery  etal.  1970) , and  the 


assumption  is  that  the  transition  from  oceanic-  to  continental- 


subduction  temperatures  on  the  fault  occurs  gradually  over  the 


period  in  which  sialic  material  begins  to  move  into  the  shear 


zone . 


In  calculating  the  heat  lost  into  the  slab  as  a result  of 


friction,  we  assume  that  conduction  of  frictional  heat  along 


the  slab  and  through  its  base  are  negligible.  Then  the  slab  is 


modelled  as  a semi-infinite  halfspace  with  a specified  thermal 


history  at  the  surface.  This  thermal  history  consists  of  a 


series  of  temperature  increments  AT^  above  the  zero-stress 


temperatures,  and  if  the  temperature  varies  linearly  between 


integration  points  then  the  flux  consumed  is: 


(z-j)  . 2Kvsine  i r 1-1 

slab  AzAk  i^1  { l A"i-1)  Vv Sine 


/z  . - z . 


/z.-z.-Az 


VsinO 


> (3) 


(Carslaw  and  Jaeger,  1959),  where  6 is  the  angle  of  slab  dip, 


AZ  is  the  depth  interval  between  integration  points  Z^,  and  tc  is 


thermal  diffusivity.  Above  the  slab  the  frictional  heat  is 


assumed  to  be  distributed  in  equilibrium  conduction  to  the  surface 


(*  t depths  less  than  the  scale  distance) , or  to  be  working  its  way 
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into  the  overriding  plate  according  to  the  model  of  a half- 
space whose  surface  temperature  is  changed  by  in  the  time 

required  to  subduct  400  km: 


Q(Zj> 

up 


ATjK 

2 


fnV~ 

Jk  (< 


k (400  km) 


(4) 


By  adding  (3)  and  (4)  we  find  the  flux  that  friction  must 
produce  to  cause  melting,  and  by  (2)  we  obtain  the  shear  stress. 
However,  this  stress  is  not  allowed  to  exceed  the  product  of 
hydrostatic  pressure  with  a coefficient  of  friction  of  0.6 
(Byerly,  1966). 

The  results  in  Figure  5 indicate  that  the  shear  stress  is 
a strong  function  of  plate  velocity  and  a weak  function  of  the 
melting  geotherm.  As  expected  from  (3),  (4)  and  (5)  the  shear 
stresses  at  1 cm/yr  are  greater  than  the  stresses  at  5 cm/yr 
by  approximately  a factor  of  the  square  root  of  the  velocity  ratio 
For  constant  K,  k,  8,  and  z,  to  maintain  melting  we  have: 


t (z) 


* *U'0> 


v-1/2 


(5) 


Stress  decreases  with  velocity,  or  conversely  that  there  may  be 
a critical  plate  velocity  below  which  continental  subduction  is 
prevented  by  friction.  The  importance  of  this  effect  depends 
critically  on  the  density  anomaly  of  the  cold  oceanic  slab  which 
provides  the  driving  force. 


— - 
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This  density  anomaly  resulted  from  the  initial  thermal 
contraction  of  the  mantle  material  forming  the  oceanic  plate 
when  it  cooled  near  a mid-ocean  ridge.  The  amount  of  contraction 
can  be  estimated  from  the  isostatic  balance  of  plates.  The 
elevation  with  respect  to  old  sea  floor  of  continents  is  5 to 
6 km  and  of  mid-ocean  ridges  is  3 to  4 km.  Noting  that  mid- 
ocean ridges  give  a reasonable  sample  of  the  asthenosphere  with 
which  we  want  to  compare  the  slab,  we  find  that  the  weight 
excess  of  a unit  section  of  oceanic  slab  is  about  8 x 108  dyne/cm2 
with  respect  to  the  asthenosphere  and  that  the  figure  for  the  net 
mass  deficiency  of  continental  lithosphere  is  about  -5  x 108  dyne/cm2. 
Unless  these  estimates  are  in  gross  error,  the  comparable  buoyancy 
of  continental  plates  and  negative  buoyancy  of  oceanic  plates 
with  respect  tc  the  asthenosphere  implies  that  hundreds  of 
kilometers  of  continental  crust  could  be  subducted  if  the  plate 
behaved  as  a unit. 

However,  the  integrity  of  the  plate  is  probably  not  preserved. 
Whatever  the  driving  force  is  in  oceanic  subduct ion,  the  devia- 
toric  stress  field  at  shallow  depths  near  the  fault  plate  must  be 
vertical  tension  and  horizontal  compression.  The  additional 
deviatoric  stress  produced  by  the  buoyancy  of  subducted  continental 
crust  would  be  vertical  tension  below  the  crust  and  vertical 
compression  above.  This  extra  stress  would  be  greatest  and  would 
interfere  constructively  with  the  driving  stresses  near  the  lower 
tip  of  the  crust.  Subsequent  shearing  would  thus  tend  to  occur 
within  the  subducting  plate  at  the  base  of  the  subducted  crust. 
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If  the  new  fault  dipped  in  the  same  direction  as  the  old 


fault,  it  would  outcrop  in  the  downing  plate  away  from  the  old 


convergence  zone-  In  the  case  of  an  island  arc-continent 


collision,  a new  fault  conjugate  to  the  old  fault  would  permit 


the  subduction  zone  to  reverse  polarity.  In  either  case  the 


slab  might  detach  ac  the  base  of  the  subducted  continental  crust 


along  the  old  fault  and  a conjugate  fault.  Hot  material 


upwelling  from  below  to  replace  the  slab  would  provide  a large  heat 


source, 


Another  means  of  approaching  this  problem  is  to  consider 
the  horizontal  convergence  of  the  two  plates  to  be  dictated  by 


distant  driving  forces  which  impose  a displacement  boundary 
condition.  Then  by  focussing  only  on  the  local  forces  which 


act  at  shallow  depths  to  oppose  the  subduction  proves*,  one  can 


seek  the  geometry  of  underthrusting  which  will  minimize  these 
forces.  This  must  be  the  geometry  of  the  fault  which  dominates 


the  subduction  process,  because  it  is  the  one  that  allows  the 


greatest  displacement  or  displacement  rate, 


The  major  forces  opposing  subduction  arise  from  thret 


effects:  friction  between  one  lithospheric  plate  and  the  other 


on  the  active  fault,  viscoelastic  resistance  to  the  bending  of 


the  downgoing  plate,  and  the  bouyancy  of  continental  crust. 


In  the  case  of  oceanic  subduction,  the  third  force  is  absent 


and  a balance  must  exist  between  the  other  two.  That  is, 


subduction  does  not  occur  through  one  plate  turning  straight 


down  because  this  would  imply  tremendous  strain  and  viscoelastic 
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resistance  in  that  plate.  Subduction  also  does  not  occur  at 
extremely  shallow  angles  because  then  the  frictional  stress 
would  act  over  a tremendous  fault  area.  At  some  intermediate 
agle  of  dip  the  total  resistance  is  minimized,  and  subduction 
proceeds  stably  over  long  intervals. 

When  continental  convergence  begins  there  are  two  changes 
which  affect  this  balance.  The  first  is  the  change  frc.u  the 
bending  of  oceanic  to  continental  lithosphere,  and  this  takes 
place  immediately.  Continental  shield  lithosphere  may  have  a 
thickness  of  about  140  km  (Toksttz  et  al . , 1967)  or  about  twice 
as  thick  as  old  oceanic  plates  (Forsyth,  1973).  Thus  it  is 
considerably  more  resistant  to  bending  and  shearing.  The  same 
bending  moment  will  produce  only  one-eighth  of  the  previous 
elastic  bending,  and  the  same  shear  will  produce  only  one-half 
the  previous  viscous  strain  rate.  This  effect  will  force  an. 
increase  in  the  radius  of  curvature  of  the  lithospheric  bend, 
so  that  downwarping  begins  in  front  of  the  old  trench  position 
and  thrusting  proceeds  at  a shallower  dip.  This  new  regime 
could  develop  soon  after  convergence  if  a new,  more  shallow- 
dipping fault  were  to  form  in  the  leading  edge  of  the  plate, 
as  discussed  above. 

The  second  change  in  the  force  balance  develops  gradually, 
as  continental  crust  on  top  of  the  subducting  plate  is  carried 
down  below  its  normal  isostatic  level.  Its  buoyancy  then  acts 
to  oppose  further  convergence  and  downward  motion,  and  the 


magnitude  of  this  buoyancy  grows  linearly  with  time.  Thus  it 
becomes  increasingly  advantageous  to  form  new  and  shallower 


thrust  faults  which  require  less  vertical  displacement  of  the 
subducted  crust  for  the  same  amount  of  horizontal  convergence. 

As  these  faults  form  the  subducted  crust  may  become  detached 
from  the  underlying  plate  and  be  itself  underthrust  as  convergence 

continues . 

We  have  shown  that  in  the  case  of  continental  plates  converging 
at  a high  velocity,  continental  subduction  at  shallow  dips  is  a 
strong  possibility.  The  resisting  forces  of  gravity,  friction, 
and  bending  could  be  balanced  by  the  known  driving  mechanism  of 
sinking  of  dense  oceanic  lithosphere.  There  is  the  interesting 
suggestion  that  convergence  of  continents  will  cause  both  the 
position  and  dip  of  the  active  inter-plate  thrust  plane  to  change. 
Much  further  work  will  be  required  before  the  exact  mechanisms  can 
be  predicted  in  proper  sequential  order,  but  the  stress,  density 
and  temperature  constraints  derived  here  will  be  essential  in 
any  attempt  to  predict  this  behavior.  As  an  intermediate  step, 
we  will  begin  consideration  of  individual  mountain  ranges  where 
continental  subduction  may  be  occurring,  and  test  the  geologic 
and  geophysical  data  available  against  the  predictions  of  models 

f.uch  as  these. 


•> 
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Tec  t o n ic^__of__the_Za2ros^ 

' The  evidence  of  marine  paleomagnetism  strongly  suggests 

that  the  Zagros  region  of  Iran  is  the  youngest  continental 
convergence  zone  on  the  Earth.  The  relative  motion  of  Arabia 
and  central  Iran  can  be  calculated  by  vector  addition  of  the 
motion  of  North  America  away  from  Europe  to  the  motion  of 
Africa  with  respect  to  North  America  (Pitman  and  Talwam,  1972) 
plus  a final  addition  of  the  opening  motion  at  the  Red  Sea. 

This  we  interpret  from  the  results  of  Phillips  (1970)  to  have 
occurred  at  2 cm/yr  for  the  last  1.5  m.y. , giving  a final  sum 
of  4.7  cm/yr  convergence  about  a pole  near  37°W,  25°N  for  the 
Zagros.  Comparison  of  this  figure  with  reasonable  estimates 
of  the  total  crustal  shortening  in  the  range  implies  a Pleistocene 

age  for  the  collision  of  Arabia  and  Iran. 

The  g.  ology  in  the  Zagros  region  is  described  by  StOcklin 
(1968)  and  Haynes  and  McQuillan  (1974).  Haynes  and  McQuillan 
present  extensive  evidence  to  show  that  the  Zagros  Crush  Zone 
which  marks  the  northeast  limit  of  deformation  is  the  site  of 
a Tertiary  oceanic  subduction  zone,  at  which  this  portion  of  the 
Tethys  Ocean  was  consumed.  The  presence  of  andesitic  volcamsm 
from  the  Eocene  to  the  Recent  along  a line  some  100  to  200  km 
Northeast  of  the  Crush  Zone  (Figvre  6)  indicates  clearly  the 
polarity  of  the  zone,  with  the  oceanic  slab  dipping  beneath 
central  Iran.  Thus,  the  main  folded  belt  of  the  Zagros  represents 
the  former  continental  shelf  of  Arabia,  with  its  thickening  wedge 
of  3 to  5 km  of  post-Jurassic  sediments  and  evaporites. 
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Many  geologic  interpretations  of  this  range  have  ascribed 
the  220  km-wide  band  of  folding  and  faulting  as  the  result  of 
decollement  on  the  layers  of  salt  as  the  basement  of  the  range 
is  thrust  under  Iran  at  the  Crush  Zone.  Chappie  (1975)  has  shown 
that  this  is  mechanically  feasible  if  the  salt  layers  are  contin- 
uous. While  this  interpretation  suggested  by  Dewey  and  J.  Bird 
(1970) , meshes  nicely  with  our  simple  theoretical  models,  it 
fails  to  explain  the  distribution  of  seismicity,  the  regional 
gravity  anomaly,  or  the  present  vertical  tectonics.  These 
factors  can  only  be  explained  if  there  are  two  major  planes  of 

underthrusting  rather  than  one. 

The  pattern  of  seismicity  as  relocated  by  Nowroozi  (1971) 
shown  in  Figure  6 is  concentrated  in  front  of  the  CrusI  Zone 
rather  than  behind  it,  where  events  would  occur  if  Arabia  were 
being  subducted  from  this  line  toward  the  Northeast.  These 
earthquakes  cannot  be  attributed  to  plate  bending  because 
Nowroozi  (1972)  and  Akasheh  (1973)  have  found  them  to  have 
nodal-plane  solutions  indicating  shallow  thrusting.  Also,  the 
pattern  of  Bouguer  gravity  anomalies  reported  by  Wilcox  et  al. 
(1972)  (Figure  7)  shows  that  the  Zagros  gravity  low  begins  from 
as  far  Southwest  as  the  edge  of  the  Persian  Gulf.  As  Figures  2 
and  3 illustrate,  it  is  impossible  to  produce  any  such  anomaly 
with  only  a slab  dipping  from  the  Crush  Zone.  Finally,  there 
is  the  existence  of  the  Mesopotamian  Trough  including  the 
Persian  Gulf  to  be  explained.  This  linear  foredeep,  which 


■ 
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thickens  towards  the  Northeast  and  is  actively  subsiding,  forms 
an  elevation  couple  with  the  Zagros  Range  that  indicates  a 

reverse  fault  on  the  NE  margin  of  the  Trough. 

Therefore,  our  preferred  model  of  Zagros  tectonics  involves 

two  major  thrust  planes  surrounding  an  internally-deformed 
median  block.  A predominantly  oceanic  slab  of  lithosphere  with 
oceanic  crust  and  sedimentary  deposits  has  been  subducted  at 
the  trench  location  now  represented  by  the  Crush  Zone  (Figure  8). 

The  oceanic  slab  is  the  former  floor  of  the  Tethys  Ocean,  and 
the  sediments  are  those  of  the  former  Arabian  continental  rise. 

The  presence  of  active  volcanoes  behind  this  line  in  Iran 

* 4-v.ie  nrpanic  slab  is  continuing  or  only 
indicates  that  motion  of  this  oceanic 

recently  halted. 

AS  the  sediments  of  the  main  range  are  of  continental  shelf 
origin,  they  overlie  a normal  continental  crust  of  perhaps  33  km 
thickness.  This  crustal  block  makes  the  Zagros  continental 
lithosphere  buoyant  as  well  as  rigid  and  it  would  therefore 
resist  subduction  along  the  steeply-dipping  plane  of  oceanic 
subduction.  instead,  this  block  appears  to  have  absorbed  about 
50  km  of  continental  convergence  by  internal  deformation  along 
numerous  minor  thrust  faults.  This  shortening,  which  would 
result  in  crustal  thickening,  is  not  expressed  at  the  surface  by 
thrusts  because  of  the  independent  mobility  of  the  sediments  above 
the  lubricating  salt  layers.  Instead,  the  competent  carbonate 
layers  form  draped  folds  over  each  major  basement  fault. 
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Primarily  to  accomodate  the  greater  flexural  rigidity  of 
the  Arabian  shield  lithosphere,  the  locus  of  tectonic  deformation 
has  moved  further  Southwest  in  the  Pleistocene.  A major  low- 
angle  thrust  dipping  NE  has  formed  along  the  former  hinge  line 
of  the  Arabian  continental  shelf,  and  this  is  responsible  for 
the  formation  of  the  Mesopotamian  Trough.  The  inferred  fault 
would  surface  beneath  a deep  layer  of  Recent  alluvium,  but  its 
location  may  be  traced  along  the  boundary  of  the  stable  area 
on  the  Tectonic  Map  of  Iran  (Stticklin  and  Nabavi,  1973)  which 
coincides  with  the  SW  limit  of  seismicity  and  the  edge  of  the 

Persian  Gulf  (Fig.  6). 

In  order  to  model  this  orogenic  zone,  information  is 
needed  on  the  dips  of  both  slabs.  This  was  obtained  from  a 
compilation  of  earthquake  locations  by  Nowroozi  (1971)  shown 
in  Figure  9.  Relatively  few  events  are  associated  with  the 
oceanic  slab,  which  may  be  inactive  or  detached,  but  they 
indicate  a dip  near  30°.  The  majority  of  seismicity  forms  a 
diffuse  band  dipping  about  15"  from  the  Persian  Gulf  margin 
where  the  latest  thrusting  is  occurring.  The  apparent  length 
of  this  band  does  not  indicate  the  amount  of  continental  sub- 
duction,  however,  because  any  new  major  thrust  is  expected  to 

be  seismic  along  its  whole  surface. 

For  lack  of  specific  information,  we  must  assume  that  the 

lithosphere  in  this  model  is  of  typical  oceanic  type.  The 
oceanic  slab  is  taken  to  be  70  km  thick;  model  parameters  are 
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given  in  Table  1.  It  could  reasonably  be  thicker  if  it  is  as 
old  as  Permian  (Dewey  et  al.  1973) . The  continental  part  of 
the  plate  is  thickened  to  140  km  (Knopoff  and  Fouda  1975)  and 
has  a low  radioactivity,  required  to  match  the  average  heat 
flow  value  of  0.88  HFU  reported  by  Coster  (1947). 

The  andesite  volcanism  (Figure  6)  behind  the  Crust  Zone 
is  caused  either  by  frictional  heating  or  by  secondary  back- 
arc  type  convection  (Andrews  and  Sleep,  1974).  As  we  are  unable 
to  model  the  latter,  we  have  set  the  frictional  stress  on  the 
subducting  oceanic  plate  at  4 kb  (the  maximum  shear  strength  of 
olivine  from  Kohlstedt  et  al.  1975)  from  40-90  km  depth  during 
the  initial  period  when  it  was  moving  at  2.7  cm/yr.  This 
creates  a high-temperature  zone  in  the  Iranian  plate  beneath 

the  volcanic  line  (Figure  8). 

In  the  second  stage  of  the  model  the  oceanic  plate  probably 
carries  traces  of  continental  rise  sediments  into  the  subduction 
zone.  From  this  point  we  have  controlled  the  stresses  so  that 
the  melting  point  of  granite  is  not  exceeded.  Melting  occurs 
below  40  km  depth,  and  this  provides  an  explanation  of  the 
synorogenic  granites  reported  by  Wells  (1969)  in  the  Crush  Zone, 
although  considerable  erosion  or  magma  upwelling  must  be  assumed 
to  reconcile  the  occurrence  and  the  source. 

About  one  million  years  ago  this  slab  stopped  cr  slowed 
considerably  and  deformation  moved  to  the  southwest.  At  this 
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point  we  introduce  a new  continental  slab  dipping  from  the 
Persian  Gulf.  Although  the  total  convergence  velocity  increased 
to  4.7  cm/yr  before  this  time,  the  new  slab  is  assumed  to  subduct 
at  only  2.6  cm/yr  because  of  simultaneous  crustal  shortening  to 
the  NE.  Melting  has  not  had  time  to  develop  on  this  new  fault, 
so  the  stress  limits  of  the  previous  section  do  not  apply. 

In  rite  ad  we  set  frictional  stress  on  the  upper  parts  of  the 
fault  equal  to  hydrostatic  pressure  multiplied  by  a coefficient 
of  friction  of  0.4  in  order  to  match  the  heat  flow  at  the 
surface.  Deeper  parts  of  the  fault  are  again  limited  to  4 kb, 
although  these  stresses  cannot  be  resolved  from  surface  heat  flow. 

This  heat  flow  was  computed  from  oil-well  temperature 
profiles  in  3-5  Y*  deep  wells  supplied  by  the  National  Iranian 
Oil  Company.  Only  smooth  mono tonic  profiles  were  used.  Heat 
flow  was  determined  by  least-squares  matching  of  these  points  to 
theoretical  temperature  profiles  calculated  using  the  conductivity 
values  and  derivatives  of  Asmari  Limestone  reported  by  Clark  (1966) 
The  results  show  a gradual  increase  of  heat  flow  from  0.90  to 
1.15  hfu  as  the  fault  is  approached  from  the  northeast  (Figure  10). 

In  a case  like  this  where  high  temperatures  have  not  yet 
developed,  the  fault  surface  separating  the  plates  would  be  a 
natural  site  for  earthquakes,  and  it  is  possible  that  the  seismic 
waves  from  these  events  will  be  useful  in  defining  the  position 
of  the  slabs.  We  have  traced  rays  in  three  dimensions  from  a 
hypothetical  earthquake  located  at  a depth  of  20  km  on  the  fault 
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in  our  first  theoretical  model.  Using  the  temperature  field 
of  Figure  2,  mantle  P-wave  velocity  anomalies  were  calculated 
using  derivatives  appropriate  for  peridotite.  In 

the  crustal  portion,  surface  velocities  were  translated 
downward.  These  anomalies  are  superimposed  on  a modified 
Jeffreys-Bullen  velocity  structure  with  a 30  k»  crust  and  a 
low  velocity  zone.  Rays  were  traced  from  the  event  at  a variety 
of  azimuths  and  angles  of  incident  (Julian,  1970,  Sleep,  1973). 
Travel-time  residuals  were  calculated  relative  to  a set  of  rays 
traced  through  an  unperturbed  model  These  residuals  are  plotted 
in  Figure  m using  the  azimuth  and  angle  of  incidence  of  each 
ray  at  150  tan  depth,  below  the  anomalous  structures.  In  the 
annulus  of  useful  rays,  all  waves  arrive  late  because  of  the 
additional  subducted  crust.  The  maximum  delay  of  1.0  to  1.3  sec 
occurs  for  waves  travelling  in  the  downdip  direction  through  the 
slab.  This  delay  is  reduced  for  large  angles  of  incidence  in  the 
downdip  direction  by  the  presence  of  the  fast  oceanic  slab  ahead 

of  the  continental  lithosphere. 

In  the  absence  of  very  accurately  located  earthquakes  in  the 

Zagros,  we  chose  to  look  at  the  travel  time  anomalies  with  the 
inverse  procedure  using  a station  located  in  the  Zagros  region 
and  earthquakes  distributed  so  as  to  give  good  azimuthal  and 
incidence  coverage.  The  WWHSS  station  at  Shiraz  was  used  for 
this  purpose.  A plot  of  observed  travel-time  residuals  at 
Shiraz  derived  from  the  I.S.C.  Bulletins  covering  the  period  of 
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1968-69  are  shown  in  Figure  11B.  This  does  not  show  any  clear 
variation  of  travel-time  residuals  across  the  strike  of  the 
range,  so  it  is  purely  negative  evidence  that  continental 
material  has  not  yet  been  subducted  beneath  Shiraz  on  a xarge 

scale. 

The  most  important  evidence  controlling  our  model  is  the 
distribution  of  Bouguer  gravity  anomalies,  because  these  princi- 
pally reflect  the  presence  of  excess  sialic  crust  below  the  33  km 
reference  thickness.  Data  compiled  by  the  U.S.  Air  Force  for 
the  Zagros  between  46°  and  56°E  is  shown  in  Figure  10  in  the 
form  of  one-degree  averages  projected  parallel  to  the  Crust  Zone. 
Some  of  the  scatter  is  due  to  the  reduction  of  the  anomaly  by 
edge  effects  near  the  ends  of  the  range,  so  our  two-dimensional 
model  will  be  fitted  to  the  maximum  anomaly  amplitudes  in  the 

central  (Shiraz)  cross-section. 

This  profile  is  contaminated  by  the  regional  -130  mg  anomaly 

of  the  Central  Iranian  plateau,  which  probably  predates  the 
collision.  To  remove  this  effect,  we  assume  that  the  plateau 
is  compensated  by  a thi;v  40  km  crust  which  extends  up  to  the 
NE  dipping  plane  of  oceanic  subduction,  and  calculate  the  edge 
effects  accordingly.  We  also  illustrate  a probable  overestimate 
of  the  gravitational  effect  of  the  Mesopotamian  Trough  (assuming 
a 600  m downwarp) . The  remainder  of  the  anomaly  must  be  due  to 
crustal  thickening  in  the  Zagros. 
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However,  the  gradients  of  this  anomaly  are  too  low  to  be 
matched  to  a dipping-slab  model  like  those  of  Figures  2 and  3. 

The  best  fitting  dip  of  the  Moho  is  actually  only  3°,  so  that 
a normal  33  km  crust  beneath  the  Persian  Gulf  passes  Northeastward 
into  a maximum  49  km  thickness  where  the  edge  of  Arabia  is 
underthrust  (Figure  8) . We  interpret  this  as  a result  of  a 
shortening  of  the  Zagros  basement  by  about  75  km  over  the  last. 

1.6  m.y . 

This  strain  is  most  intense  on  the  NE  margin  where  the 
crustal  thickening  and  surface  deformation  are  greatest. 

According  to  the  fault-plane  solutions  the  present  mechanism  is 
shallow  underthrusting,  so  the  Moho  is  moving  NE  more  rapidly 
than  down.  A slip  of  27  km  (26  km  NE  and  7 km  down)  on  the  new 
boundary  thrust  would  be  consistent  with  the  seismic  and  gravity 
data  in  the  SW  half  of  the  Zagros. 

This  model  has  been  derived  under  the  assumption  of  plate 
tectonics,  which  provides  the  necessary  limitations  on  the  number 
of  variable  parameters.  Of  course,  if  continental  lithosphere 
does  not  have  the  strength  to  deform  in  large  coherent  units,  then 
any  number  of  non-unique  explanations  can  be  attached  to  the 
available  geologic,  gravity,  and  seismicity  data.  We  are  supported 
in  our  assumption  by  the  similarity  of  the  Zagros  structures  over 
the  great  length  of  the  range  which  have  evolved  to  a common 
pattern  despite  probable  initial  variations  in  coastline  shape 
and  crustal  thickness.  The  smooth  downwarp  of  the  Arabian 


lithosphere  over  150  to  400  kin  as  it  approaches  the  Mesopotamian 
Trough  is  additional  evidence  of  lithospheric  strength  and  the 
applicability  of  plate  tectonic  concepts. 

This  model  implies  two  important  conclusions  about  the 
continental-convergence  process.  First,  that  subduction  r.iy 
not  occur  at  the  old  trench  site  as  in  our  simple  models  and 
those  of  Dewey  and  Bird.  Instead  the  crust  begins  to  deform 
and  thicken.  Second,  that  the  viscoelastic  resistance  to  the 
bending  of  the  thick  lithosphere  of  continental  shields  may  be 
sufficient  to  force  the  formation  of  a new  fault  at  a shallower 
angle.  As  subduction  begins  along  this  new  fault  it  will 
isolate  the  continental  shelf  block  from  further  subduction  and 
result  in  its  preservation  at  the  core  of  a future  mountain  range. 
Thus  the  Zagros  Mountains  may  hold  the  key  to  the  final  reconcili- 
ation of  the  classical  and  plate-tectonic  definitions  of  the 
geosyncline . 

Conclusions 

Theoretical  models  calculated  in  this  paper  strongly  depend 
on  input  parameters  which  cannot  be  specified  uniquely.  For  a 
reasonable  set  of  parameters,  it  can  be  stated  that  during  the 
subduction  process  there  is  insufficient  time  for  radioactive 
heating  or  heat  conduction  through  the  bottom  of  the  lithot  phere 
to  significantly  affect  crustal  temperatures.  These  processes 
become  important  after  subduction  has  stopped.  Calculations  show 
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that  after  30  m.y.  of  quiescence  a continental  slab  of  150  km 
length  and  2.5  x 10"6  erg/gm-sec  average  radioactive  heat 
production  will  be  raised  to  high-grade  metamorphic  temperatures. 
After  an  additional  70  m.y.  melting  would  develop  in  the  base  of 
the  slab,  leading  to  plutonism.  The  rate  of  this  regional- 
metamorphic  evolution  is  directly  controlled  by  crustal 
radioactivity  and  surface  erosion  rates,  so  that  there  is  wide 
variability  possible  in  the  timing  and  severity  of  the  resulting 
orogeny.  However,  the  rapidity  of  observed  metamorphic  events 
suggests  an  influx  of  heat  from  the  mantle.  Frictional  heating 
could  be  an  important  mechanism  in  thermal  regime  of  the  shallow 
part  of  the  subduction  zone. 

The  highest  temperature  reached  through  frictional  heating 
is  controlled  by  the  profile  of  shear  stress  versus  depth  and  the 
square  root  of  the  plate  convergence  velocity.  At  a velocity  of 
1 cm/yr  and  a shear  stress  of  1 kb  frictional  heating  is  too  weak 
to  cause  any  major  change  in  the  low-temperature  anomaly  introduced 
by  the  slab.  However,  at  6 cm/yr  with  similar  stresses  a high- 
temperature  zone  is  formed  on  top  of  the  subducting  crust,  and 
temperatures  may  reach  the  solidus  of  granite.  This  will  result 
in  minor  igneous  activity,  but  is  primarily  important  because  it 
sets  a limit  on  the  frictional  stress  at  each  dept'i. 

At  5 cm/yr,  these  stress  limits  have  been  found  to  fall 
below  1 kilobar  at  depths  greater  than  25  km.  Stick-slip  behavior 
passes  into  melting  or  creep  strain  below  about  15  km  depth. 
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Thus  the  seismicity  of  continental  convergence  zones  rs 
necessarily  shallow  except  for  events  within  the  slab.  At 
lower  velocities  the  stresses  required  to  produce  melting  are 

unreasonable . 

The  other  important  forces  in  continental  convergence 
result  from  crustal  buoyancy  and  resistance  to  the  bending  of 
the  lithosphere.  Both  tend  to  force  the  angle  of  subductron 
to  become  more  shallow  with  time.  While  the  effect  of  buoyancy 
builds  up  gradually  with  the  amount  of  crust  subducted,  the 
greater  thickness  of  continental  shield  lithosphere  changes  the 
mechanism  of  plate  bending  as  soon  as  the  continents  meet. 
Subduction  on  the  old  oceanic  plane  is  strongly  resisted,  and 
the  formation  of  shallower  thrust  faults  is  required. 

This  process  is  shown  to  be  active  in  the  Zagros  Mts.,  where 
there  is  a very  young  continental  convergence  zone.  There 
mountains  are  developing  on  a horizontal  triangular  prism  of 
detached  crust  between  the  old  oceanic  subduction  plane  to  the 
HE  and  a new  shallower  fault  plane  to  the  SW.  Continental 
underthrusting  is  continuing  at  a different  site  from  the 
continental  suture  line.  The  crust  of  the  Zagros  range,  which 
is  the  old  Arabian  continental  shelf,  is  being  thickened  by 
underthrusting  in  the  first  stage  of  the  erogenic  process. 


, 


ACKNOWLEDGEMENTS 


I 


We  thank  John  Minear  for  his  work  in  the  original 
development  of  the  temperature  program.  Seth  Stein  developed 
and  assisted  with  programs  used  in  tracing  seismic  rays. 

The  National  Iranian  Oil  Company  courteously  provided 
temperature  data  from  eighteen  of  their  oilfields.  The 
Aeronautical  Chart  and  Information  Center  of  the  U.S.  Air 
Force  provided  a compilation  of  available  gravity  data. 

This  research  was  supported  by  the  Advanced  Research 
Projects  Agency  and  monitored  by  the  Air  Force  Office  of 
Scientific  Research  under  contract  F44620-71-C-0049  and 
NATO  Science  Grant  #568  with  Massachusetts  Institute  of 
Technology  and  by  Grant  DEF  74-22338  from  the  National 
Science  Foundation  to  Northwestern  University. 


-43- 


REFERENCES 

Akascheh,  B.  , Mechanism  of  the  earthquake  of  April  10,  1972  in 

Qir  (Iran),  Zeitschrift  fur  Geophysik,  39 , 1055-1061,  1973. 

Andrews,  D. J. , and  N.H.  Sleep,  Numerical  modelling  of  tectonic 

flow  behind  island  arcs,  Geophys.  J.  Roy.  Astron^ §°.c_l> 

237-251,  1974. 

Brown,  G.C. , andW.S.  Fyfe,  The  production  of  granitic  melts 
during  ultrametamorphism,  Contr.  Mineral.  Petrol.,  20, 
310-318,  1970. 

Byerly,  J.D.,  The  frictional  characteristics  of  Westerly 

granite,  Ph.D.  Thesis,  Massachusetts  Institute  of  Technology, 
Cambridge,  Mass.,  1966. 

Carslaw,  H.S.,  and  J.C.  Jaeger,  Conduction  of  Heat  in  Solids, 
Clarendon  Press,  Oxford,  pp.  62-75,  1959. 

Chappie,  W.M.,  Mechanics  of  thin-skinned  fold  and  thrust  belts, 
Trans.  Am.  Geophys.  Uni^n,  56 , 457,  1975. 

Clark,  S.P.,  and  A. E.  Ringwood,  Density  distribution  and 

constitution  of  the  mantle,  Rev.  Geophys.,  2_,  35-88,  1964. 

Clark,  S.P.,  Jr.,  Thermal  conductivity  in  Handbook  of  Physical 
Constants,  Geol.  Soc.  Amer.  Mem.,  97,  459-432,  1966. 

Coster,  H.P.,  Terrestrial  heat  flov*  in  Persia,  Monthly  Notices, 
Geophysical  Supplement  5,  R.  astr . Soc . , 131-145,  1947. 

Dana,  J.D.,  On  some  results  of  the  earth's  contraction  including 
a discussion  of  the  origin  of  mountains  and  the  nature  of 
the  earth's  interior,  Amer.  J . Sci. , 5_/  No,  5,  423-443, 
474-475;  No.  6,  6-14,  104-115,  161-172,  304,  381-383,  1873. 


— - ~ r ' ~ - ***- - 


Y 


-44- 


J 


Dewey,  J.F.,  and  J.M.  Bird,  Mountain  belts  and  the  new  global 
tectonics,  J.  Geophys.  Res.,  75,  2625-2647,  1970. 

Dewey,  J.F. , W.C.  Pitman,  III,  W.B.F.  Ryan,  and  J.  Bonnin, 

Plate  tectonics  and  the  evolution  of  the  Alpine  system, 

Geol.  Soc.  Amer.  Bull.,  84 , 3137-3180,  1973. 

Emery,  K.O.,  E.  Uchupi,  J.D.  Phillips,  C.O.  Bowin,  E.T. 

Bunce,  and  S.T.  Knott,  Continental  rise  off  eastern  North 
America,  Amer.  Assoc.  Petrol.  Geol.  Bull.,  5_4,  44-108,  1970. 

Fitch,  T.J.,  Earthquake  mechanism  in  the  Himalayan,  Burmese, 

Andaman  regions  and  continental  tectonics  in  Central  Asia, 

J.  Geophys.  Res.,  75,  2699-2708,  1970. 

Forsyth,  D.W.,  Anisotropy  and  structural  evolution  of  the  oceanic 
upper  mantle,  Ph.D.  Thesis,  Massachusetts  Institute  of 
Technology  and  Woods  Hole  Oceanographic  Institution, 
Cambridge,  Mass.,  1973. 

Hall,  James,  On  the  formation  of  mountains,  Canadian _J^,  5, 
542-544,  1860. 

Haynes,  S.J.,  and  H.  McQuillan,  Evolution  of  the  Zagros  suture 
zone,  southern  Iran,  Geol.  Soc.  Amer.  Bullw  85,  739-744, 

1974. 

Julian,  B.,  Ray  tracing  in  arbitrarily  heterogeneous  media, 

Technical  Note  1970-45,  Lincoln  Laboratory,  Massachusetts 
Institute  of  Technology,  Lexington,  Mass.,  1970. 

Knopoff,  L.  and  A. A.  Fouda,  Upper-mantle  structure  under  the 
Arabian  Peninsula,  Tectonophysics , 26,  121-134,  1975. 


Kohlstedt,  D.L.,  C.  Goetze,  and  W.B.  Durham,  Experimental 

deformation  of  single  crystal  olivine  with  application  to 
fl jw  in  the  mantle,  preprint,  1975. 

LePichon,  X.L.,  Sea-floor  spreading  and  continental  drift, 

J.  Geophys.  Res.,  7 3 , 3661-3697,  1968. 

MacDonald,  G.,  Calculations  on  the  thermal  history  of  the  earth, 

J.  Geophys.  Res.,  64 , 1967-2000,  1959. 

Minear,  J.W.,  and  M.N.  ToksOz,  Thermal  regime  of  a downgoing 

slab  and  new  global  tectonics,  J.  Geophys.  Res.  , 75.,  1397- 
1419,  1970a. 

Minear,  J.W.  and  M.N.  ToksOz,  Thermal  regime  of  a downgoing 
slab,  Tectonophysics , 10 , 367-390,  1970b. 

Minear,  J.W.,  and  M.N.  ToksOz,  Reply,  J.  Geophys.  Res.  , 7j6, 
610-612,  1971. 

Naylor,  R.S.,  Acadian  orogeny:  an  abrupt  and  brief  event, 

Science , 172 , 558-559,  1971. 

Nowroozi,  A. A. , Seismo-tectonics  of  the  Persian  Plateau, 

Eastern  Turkey,  Caucasus,  and  Hindu-kush  regions,  Bull^ 
Seismo.  Soc.  Amer . , 61,  317-341,  1971. 

Nowroozi,  A. A.,  Focal  mechanism  of  earthquakes  in  Persia, 

Turkey,  West  Pakistan,  and  Afghanistan  and  plate  tectonics 
of  the  Middle  East,  Bull.  Seismo.  Soc.  Amer.,  62,  823-850, 

1972. 

Oxburgh,  E.R.,  and  D.L.  Turcotte,  Thermal  gradients  and  regional 
metamorphism  in  overthrust  terrains  with  special  reference 
to  the  Eastern  Alps,  Schweiz.  Mir.  Petr.  Mitt.,  5_4,  641-662, 


1974. 


Phillips,  J.D.,  Magnetic  anomalies  in  the  Red  Sea,  Phil^ 

Trans.  Roy.  Soc.  Lon- * 267 .»  205-217,  1970. 

Pitman,  W.C.,  and  M.  Talwani,  Sea-floor  spreading  in  the 

North  Atlantic,  Geol.  Soc.  Amer.  Bullw  8^,  617-646,  1972. 

Quereshy,  M.N.,  S.  Venkatachalam,  and  C.  Subrahmanyam,  Vertical 
tectonics  in  the  Middle  Himalayas:  an  appraisal  from  recent 
gravity  data,  Geol.  Soc.  Amer.  Bull..,  85,  921-926,  1974. 

Skinner,  B.J.,  Thermal  expansion,  in  Handbook,^ 

Geol.  Soc.  Amer.  Mem.,  97,  75-96,  1966. 


Sleep,  N.H.,  Teleseismic  P-wave  transmission  through  slabs, 

Bull.  Seismol.  Soc.  Amer. , 63,  1349-1373,  1973. 

Sleep,  N.H.,  Stress  and  flow  beneath  island  arcs,  Geophys.  J^_ 
Roy,  astr.  Soc.,  43,  in  press,  1975. 

Smith,  A.T.,  and  M.N . ToksOz,  Stress  distribution  beneath  island 
arcs,  Geophys.  J.  Roy,  astr.  Soc._,  29,  289-318,  1972. 

Smith,  A.T.,  Time-  dependent  strain  accumulation  and  release  at 
island  arcs,  EOS  Trans.  AGU,  55,  427,  1974. 

Stftcklin,  J.,  Review  of  Iranian  tectonics,  Bull.  Amer.  Assoc^ 


Petrol.  Geol. , 52 , 1229-1258,  1968. 

Stttcklin,  J.,  and  M.H.  Navabi,  Tectonic  map  of  Iran,  Geological 

Survey  of  Iran,  Tehran,  1973. 

Toksflz,  M.N,,  M. A.  Chinnery,  and  D.L.  Anderson,  Inhomogeneities 
in  the  earth's  mantle,  Geophys.  J.  Roy-  astr.  Soc.,  13, 
31-59,  1967. 


— 1— . .. 


1 — 


Toksttz,  M.N. , 


J.W.  Minear,  and  B.  Julian,  Temperature  field 
and  geophysical  effects  of  a downgoing  slab,  J.  Geophys._ 

Res. , 76,  1113-1138,  1971. 

ToksOz,  M. N . , N.H.  Sleep,  and  A.T.  Smith,  Evolution  of  the 
downgoing  lithospnere  and  the  mechanisms  of  deep  focus 
earthquakes,  Geophys.  J.  Roy,  astr.  Soc.,  3J3,  285-310,  1973. 

Tuttle,  O.F.,  and  N.L.  Bowen,  Origin  of  granite  in  the  light  of 

experimental  studies  in  the  system  NaAlSi308  - KAlSi308-Si02-H20 
Geol.  Soc.  Amer.  Mem.  7 4 , 195  8. 

Turcotte,  D.L. , and  E.R.  Oxburgh,  A fluid  theory  for  the  deep 

structure  of  dip-slip  fault  zones,  Phys.  Earth  Planet.  Int^, 

1,  381-386,  1969. 

Walcott,  R.I.,  Isostatic  response  to  crustal  loading  in  Canada, 

Can.  J.  Earth  Sci.,  _1,  716-727,  1970a. 

Walcott,  R.I.,  Flexural  rigidity,  thickness,  and  viscosity  of 
the  lithosphere,  J.  Geophys.  Res.,  7J5,  3941-3954,  1970b. 

Watts,  A.B.,  and  M.  Talwani,  Gravity  anomalies  seaward  of  deep-sea 
trenches  and  their  tectonic  implications,  Geophys.  J.  Roy^_ 
astr.  Soc. , 36 , 57-90,  1974. 

Wells,  A. J. , The  Zagros  crush  zone  and  its  tectonic  implications, 
Geol.  Mag.,  106,  385-394,  1969. 

Wilcox,  L.E.,  W.J.  Rothermel,  J.T.  Voss,  A geophysical  geoid  of 
Eurasia,  Trans.  Am.  Geophys.  Union,  5J3,  343,  1972. 


-48 


-49- 


Figure  Captions 


Figure  1.  Three  stages  in  the  evolution  of  a hypothetical 


continental  convergence  zone:  A)  Subduction  of  the  inter- 


vening oceanic  plate.  An  Atlantic-type  margin  approaches 


a Pacific-type  trench  with  andesitic  volcanism  produced  by 


underthrusting.  B)  Subduction  of  continental  rise  sediments 


and  then  crust  of  the  Continental  shelf,  perhaps  with  granite 


magma  produced  by  frictional  heat.  C)  Formation  of  a new 


thrust  within  continental  crust  and  suoduction  beneath  the 


continental  shelf  with  further  granite  formation. 


Figure  2.|  Temperature,  heat  flow,  and  Bouguer  gravity  sections 


after  subduction  of  continental  crust  for  8 m.y.  at  30°  dip 


and  1 cm/yr.  Mafic  portion  of  the  lithosphere  below  the 


Moho  is  shaded.  Frictional  stress  on  the  fault  is  1 kb. 


Figure  3.  Temperature,  heat  flow,  and  Bouguer  gravity  sections 


after  subduction  of  continental  crust  for  2.5  m.y  at  15  dip 


and  6 cm/yr.  Frictional  stress  decreases  from  1.5  kb  at 


10  km  depth  to  1 kb  at  80  km. 


Figure  4.  Evolution  of  a convergence  zone  without  motion: 


A)  Temperature  and  heat  flow  section  evolved  fr.'m  Figure  3 


after  30  m.y.  Dashed  lines  show  former  outlines  of  the  slab. 


B)  Temperature  and  heat  flow  section  evolved  from  4 A after 


70  m.y.  additional  time.  Hatched  area  is  partially  melted 


and  a source  of  granite  magma. 
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Figure  5.  Two  different  solidus  curves  (center),  for  granite 
saturated  with  water  and  dry  granite  with  muscovite. 
Temperatures  on  fault  surface  (bottom)  with  subduction 
at  30°  dip  and  5 or  1 cm/yr  and  no  friction.  Stresses  at 
each  depth  (top)  necessary  to  raise  these  temperatures 
to  either  solidus.  Near  surface,  stress  is  limited  by 
pressure  and  the  coefficient  of  friction. 

Figure  6.  Outline  map  of  the  Zagros  after  Stttcklin  and  Nabavi 
(1973):  1.  Zagros  Crush  Zone  and  NE  limit  of  deformation; 

2.  SW  limit  of  deformation;  3.  Recent  volcanoes;  4.  Earthquake 
epicenters  1950-1565  after  Nowroozi  (1971);  5.  Tertiary 
volcanic  rocks,  mainly  andesite;  6.  Persian  Gulf  and  various 
seas  (inset) . 

Figure  7.  Gouguer  gravity  anomaly  in  the  Zagros  region  after 
U.S.A.F.  "Bouguer  gravity  anomaly  map  of  Asia"  (Wilcox 
et  al. , 1972).  Contour  interval  50  milligalls. 

Figure  8.  Preferred  model  of  present  Zagros  temperatures  (top) 
and  tectonics  (bottom)  . Computation  of  thermal  model 
described  in  text.  Diagram  shows  major  Recent  thrust,  minor 
thrusts  under  surface  folding,  subducted  continental  rise 
sediments  (dashed) , and  thickening  of  the  lithosphere  toward 
the  Arabian  shield. 

Figure  9.  Projection  of  earthquakes  relocated  by  Nowroozi  (1971) 
onto  plane  of  Figure  8 with  same  distance  scale.  Principal 
stresses  by  nodal-plane  solutions  of  Nowroozi  (1972)  and 
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Akascheh  (1973).  Dashed  lines  show  inferred  planes  of 
oceanic  (left)  and  continental  (right)  subduction. 

Figure  10.  Heat  flow  curve  from  model  of  Figure  8 compared 

with  one  heat  flow  determination  (large  dot)  and  heat  flow 
estimates  from  oil-well  temperatures  (small  dots) . At 
bottom,  total  Bouguer  gravity  anomaly  of  the  model  is 
subdivided  into  constituent  parts  and  compared  with  one- 
degree  averages  projected  onto  the  cross-section  (dots). 

Figure  11A.  Anomalous  delay  in  seconds  of  P waves  from  an 
event  at  20  km  depth  on  the  fault  in  Figure  2.  Delays 
are  plotted  according  to  ray  azimuth  and  incidence  angle 

at  a depth  of  150  km. 

11B.  Travel-time  residuals  of  less  than  3 sec  magnitude 
recorded  at  Shiraz  (Figure  8)  during  1/68-4/69  according 
to  the  Bulletin  of  the  I.S.C.  Average  residual  value  of 
-0.5  sec  has  been  removed.  Dashed  line  shows  the  strike 


of  the  range. 
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1 , 4 North  American-Eurasian  Plate  Boundary  in  Northeast_As ta 

by  M.E.  Chapman  and  S.C.  Solomon 
ABSTRACT 

The  intracontinental  portion  of  the  boundary  between  the 
North  American  and  Eurasian  plates  can  be  identified  on  the 
basis  of  seismicity,  recent  tectonics  and  earthquake  focal 
mechanisms.  The  simplest  plate  geometry  that  can  explain 
these  data  involves  a North  American-Eurasian  boundary  that 
extends  from  the  Nansen  ridge  through  a broad  zone  of  defor- 
mation in  northeast  Asia  to  the  Sea  of  Okhotsk  and  thence 
southward  through  Sakhalin  and  Hokkaido  to  a triple  junction 
in  the  Kuril-Japan  trench.  Such  a configuration  can  account 
quantitatively  for  the  slip  vectors  derived  from  earthquake 
mechanisms  in  Sakhalin  and  Hokkaido,  On  the  basis  of  new 
slip  vector  data,  the  North  American-Eurasian  angular  velocity 
vector  is  revised  only  slightly  from  previous  determinations. 
The  intracontinental  plate  boundary  is  diffuse  and  may  be 
controlled  by  ancient  plate  sutures.  Deformation  within  about 
10°  of  the  rotation  pole,  which  lies  very  near  the  boundary, 
cannot  be  modeled  by  rigid  plate  tectonics.  These  character- 
istics of  intracontinental  plate  boundaries  are  related  to  the 
greater  thickness  and  heterogeneity  of  continental  lithosphere 
and  to  the  influence  of  continents  on  the  plate  tectonic 
driving  forces. 


INTRODUCTION 


The  theory  of  plate  tectonics  is  based  on  the  idea  that  the  earth's 
surface  may  be  subdivided  into  a small  number  of  rigid  plates  (Morgan , 1968; 

Le  Pichon,  1968) . In  oceanic  areas,  the  boundaries  between  major  plates  are 
readily  defined  by  the  distribution  of  earthquakes  and  characteristic  bathymetric 
features  and  are  typically  no  more  than  a few  kilometers  in  width.  Where  plate 
boundaries  bisect  continental  masses,  however,  the  place  where  one  plate  ends 
and  another  begins  is  generally  much  more  difficult  to  locate  with  confidence. 

The  problem  of  identifying  a plate  boundary  within  a continent  is  heightened 
when  the  relative  velocity  of  the  two  plates  is  small. 

In  this  paper,  we  examine  in  detail  several  aspects  of  one  such  intra- 
continental plate  boundary:  the  boundary  in  northeast  Asia  between  the  Eurasian 
and  the  North  \merican  plates.  There  are  several  reasons  for  such  a study.  If 
the  concept  of  distinct  plates  is  valid,  then  each  plate  must  be  encircled  by 
some  closed  curve  on  the  earth's  surface.  Complete  specification  of  the  boundary 
of  each  plate  is  formally  necessary  to  conduct  certain  tests  of  driving  force 


models  for  plate  tectonics  (Solomon  rind  Sleep,  1974;  Forsyth  and  Uyeda,  1975; 
Solomon  et  al.,  1975)  . More  inter.’  -ting  are  the  underlying  causes  for 
the  diffuse,  rather  ill -defined  nature  of  intracontinental  plate  boundaries, 
including  the  possibilities  that  continental  lithosphere  is  thicker,  more 
difficult  to  push  or  pull  across  the  earth's  surface,  and  more  heterogeneous 
than  oceanic  lithosphere. 

We  begin  with  a review  of  historical  seismicity  and  recent  tectonic 
activity  in  northeast  Asia.  An  evaluation  of  several  possible  plate  config- 
urations to  explain  these  data  is  then  made.  The  simplest  explanation  compatible 
with  seismic  and  tectonic  evidence  is  that  the  present  North  American-Eurasian 
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plate  boundary  extends  from  the  Nansen  ridge  in  the  Arctic  ocean  :hrough  a 
broad,  tectonically  a .'-u  belt  in  northeast  U.S.S.R.  (Demenitskaya  and 
Karasik,  1969;  Grachev  et  al.,  1970;  Churkin,  1972)  to  the  Sea  of  Okhotsk, 
and  thence  southward  through  Sakhalin  and  Hokkaido  to  a triple  junction  at  the 
Japan-Kurile  trench.  Earthquake  source  mechanisms  are  consistent  with  this 
view  except  in  the  immediate  vicinity  of  the  relative  rotation  pole  for  the 
two  plates.  Finally,  some  thoughts  are  offered  on  the  factors  affecting  the 
location  and  nature  of  intracontinental  plate  boundaries. 


SEISMICITY 

The  boundaries  between  plates  are  seismically  active.  Seismicity  is  the 
primary  basis  for  identifying  the  North  American-Eurasian  plate  boundary  in  the 
Atlantic  and  Arctic  Oceans  as  the  mid-Atlantic  and  Nansen  ridge  systems, 
respectively.  The  correctness  of  such  an  identification  is  apparent  by 
inspection  from  a global  seismicity  map  based  on  only  a few  years  of  data, 
such  as  that  in  Figure  1.  How  the  plate  boundary  continues  from  the  Arctic 
ocean  onto  the  Eurasian  continent  is  by  no  means  clear  from  the  figure,  however. 

As  an  aid  in  better  defining  this  plate  boundary  in  northeast  Asia,  we  have 
plotted  in  Figure  2 the  epicenters  of  all  instrumentally  located  shallow  seismic 
events  in  the  region  for  the  period  1909-1973.  The  original  epicenters  were  obtained 
from  a number  of  earthquake  data  sources  (Linden,  1961;  Savarensky  et  al^,  1962; 
Hodgson  et  al.,  1965;  Sykes,  1965;  Solov'yev,  1965;  Int.  Seismol.  Cent.,  1966- 
1973;  Nat.  Ocean.  Atmos.  Admin.,  1970-1973;  Int.  Geod.  Geophys.  Ihuon,  1918- 
1963;  Academy  of  Sciences  USSR,  1963-1970).  Where  possible,  the  events  prior  to 

1952  were  relocated  using  modern  travel  time  tables  (Herrin,  1968)  to  obtain  more 
accurate  locations.  Observed  P-wave  and  S-wave  arrival  times  are  from  the 
bulletin  of  the  Int.  Geod.  Geophys.  Union  (1918-1951) ; the  P-wave  and  S-wave 


travel  time  uncertainties  were  assumed  to  be  1.5  and  3.0  seconds,  respectively. 

The  relocated  events  are  listed  in  Table  1. 

Several  conclusions  may  be  drawn  from  these  seismicity  maps.  In  oceanic 
regions,  the  plate  boundary  defined  by  the  seismic  belt  is  very  narrow,  possibly 
as  little  as  10  km  in  width  (Figure  1).  However,  at  the  Eurasian  continental 
margin,  the  seismic  belt  widens  to  300  km  (Figure  2) . A broad  seismic  zone 
extends  from  the  Laptev  sea  across  northeast  U.S.S.R.  to  the  northern  Sea  of 
Okhotsk.  The  zone  is  600  km  wide  at  its  maximum  width,  and  includes  both 
large  and  small  earthquakes  throughout  its  entire  extent. 

This  broad  seismically  active  band  appears  to  terminate  abruptly  at  the 
Sea  of  Okhotsk.  The  apparent  aseismicity  of  the  Sea  of  Okhotsk  may  be  an  artifact 
of  a sparse  instrumental  network  or  may  possibly  be  real.  Until 
1968  there  were  few  seismometers  in  the  region;  the  nearest  stations  were 
Magadan  (59°  33’N,  150°  48’E),  Okha  (53°  33%  142°  56’E),  and  Yakutsk  (62°  01% 
129°  43'E).  A major  earthquake  in  the  region,  with  body  wave  magnitude  5 1/4, 
occurred onMay  10,  1947  at  57.9°N,  141. 9#E.  Prior  to  1968  there  were  no  small 
events  recorded  in  this  area.  However,  after  five  new  seismic  stations  were 
installed  in  Yakutia  in  1968,  small  events  were  detected  in  the  Sea  of  Okhotsk 
close  to  the  stations  (Academy  of  Sciences  USSR,  1968).  Because  of  this  increase 
in  reported  events  and  the  one  major  event,  it  is  likely  that  the  Sea  of  Okhotsk 
region  is  seismically  active,  even  though  earthquakes  with  magnitudes  greater 
than  4 are  rare. 

Another  feature  of  Figure  2 is  a zone  of  small  magnitude  earthquakes  from 
120°  to  approximately  135°E  longitude  at  about  56°N  latitude.  This  zone  is  an 
extension  of  the  Lake  Baikal  seismic  belt  (Gutenberg  and  Richter,  1949).  It  is 
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not  clear  from  the  seismicity  whether  this  zone  continues  eastward  to  the 

» 

eeismically  active  portions  of  the  Sea  of  Okhotsk. 

There  are  many  shallow  events  on  Sakhalin,  an  area  known  for  its  relatively 

high  seismicity  (Solov'yev,  1965) . The  north-trending  seismic  belt  on  Sakhalin 
is  about  300  km  wide,  and  appears  to  continue  through  Hokkaido  to  the  Japan- 

Kurile  trench. 

There  are  numerous  earthquakes  located  in  the  Kuril  and  Aleutian  trenches, 
and  many  small,  shallow  events  extend  inland  on  Kamchatka.  Also  on  Kamchatka 
there  appears  to  be  a short  seismic  belt  extending  northwards  from  the  Koman- 
dorski islands.  This  belt  is  probably  related  to  under thrusting  of  the  Pacific 

plate  (Cormier,  1975) . 

f I 
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RECENT  TECTONICS 


A brief  review  of  field  geological  and  geophysical  evidence  for  recent 
tectonic  activity  in  the  Sea  of  Okhotsk  region  is  a necessary  preliminary  to 
the  discussion  of  possible  plate  boundaries  in  the  area  and  to  the  interpretation 

of  earthquake  focal  mechanisms. 

The  principal  Cenozoic  tectonic  features  of  central  and  southern  Kamchatka 
are  illustrated  in  Figure  3 (Alverson  et  al. , 1967) . Most  faults  in  Kamchatka 
are  parallel  to  the  Kuril  trench,  and  where  presently  active  are  likely 
related  to  subduction  of  the  Pacific  plate  rather  than  the  plate  boundary  in 
question.  The  only  major  fault  system  that  trends  in  an  east-west  direction 
is  the  Kronoki-Krutogorova  fault  zone  (Suprenko  and  Dekin,  1968;  Suprenko  et  al.. 
1973).  In  the  Kronoki  peninsula  the  faults  in  this  zone  are  right  lateral,  with 
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a maximum  horizontal  offset  of  20-25  km.  In  the  western  portion  of  the  Kronoki- 
Krutogorova  fault  zone  the  sense  of  motion  is  left-lateral,  opposite  to  the 
sense  of  motion  in  eastern  Kamchatka.  These  western  faults  have  offsets  similar 
in  amplitude  to  those  of  the  faults  in  the  Kronoki  peninsula.  The  fault  zone 
may  form  the  southern  boundary  of  a rift  system,  active  in  Plio-Pleistocene 
times,  in  the  central  Kamchatka  basin  (Suprenko  et  al.,  1973) . 

Sakhalin  is  dominated  structurally  by  compressive  features  such  as  faults 
and  folds  that  trend  north-south  along  the  longitudinal  axis  of  the  island. 

One  of  the  primary  faults  is  the  central  Sakhalin  fault,  a thrust  fault  with 
a meridional  trend  and  a westerly  dip  of  approximately  70° . The  strike  of 
drag  folds  and  second  order  faults  indicates  some  right-lateral  movement  along 
the  main  fault.  This  fault  is  dated  as  being  active  in  the  Late  Miocene  and 
Pliocene,  and  is  still  active  today  (Zanyukov,  1971).  Quaternary  displacements 
have  been  measured  and  the  epicenters  of  crustal  earthquakes  are  located  in 
the  fault  zone. 

In  the  Schmidt  peninsula  (northern  tip  of  Sakhalin)  there  are  also  north 
to  no.thwest  trending  thrust  faults  which  exhibit  some  right-lateral  motion. 

The  sense  of  horizontal  displacement  is  indicated  by  the  inclination  of 
slickensides,  displaced  features  and  drag  folds.  There  is  up  to  14  km  of  right- 
lateral  offset.  The  main  tectonic  activity  was  in  Plio-Pleistocene  time. 
However,  the  presence  of  tectonic  scarps,  rockfalls  and  slides,  and  the 
reworking  of  stream  drainage  patterns  indicate  that  the  movements  are  continuing 
today  (Rozhdestvenskiy , 1973)  . 

That  compressive  forces  have  acted  upon  Sakhalin  is  also  indicated  by 
folding.  The  axes  of  anticlines  and  synclines  trend  north-south  in  Sakhalin 


parallel  to  the  strike  of  the  thrust  faults  (Pushcharovskiy , 1965;  Gal  tsev^ 
Bezyuk,  1968).  During  the  latest  episode  of  folding,  denoted  the  Sakhalin 
Orogeny,  Pliocene  deposits  were  folded.  This  episode  of  folding  has  continued 
until  the  present.  Geodetic  measurements  indicate  that  vertical  crustal  move- 
ments of  3 to  9 ram/year  are  occuring  in  southern  Sakhalin  (Zakharov  and 
Yakushko,  1972). 

II 

SOME  POSSIBLE  PLATE  CONFIGURATIONS 

There  are  a number  of  ways  northeast  Asia  may  be  subdivided  into  plates 
consistent  with  the  seismicity  and  recent  tectonic  activity  discussed  above. 

All  of  these  plate  descriptions  are  to  some  extent  inadequate  for  they  fail 
x to  account  fully  for  intraplate  deformation  and  the  finite  width  of  intra- 

continental plate  boundaries.  Nonetheless,  we  prefer  one  such  description 
that  is  at  the  same  time  simple,  in  approximate  agreement  with  the  evidence 
outlined  in  preceding  sections,  and  useful  for  making  quantitative  predictions. 

The  primary  test  of  a proposed  plate  model  is  whether  the  sense  of 
. ! motion  at  plate  boundaries  predicted  by  the  relative  angular  velocity  vector 

of  the  adjacent  plates  is  consistent  with  seismic  and  tectonic  evidence.  To 
make  this  test  we  need  reasonably  accurate  estimates  of  the  Eurasian-North 
American  rotation  pole  and  rate.  A number  of  determinations  of  theee 
quantities  have  been  made  by  a variety  of  techniques;  these  are  summarized  in 
Table  2.  Probably  any  of  the  three  most  recently  published  solutions  are 
suitable  for  the  purposes  of  this  section.  In  a later  section  we  derive  a new 

* 

Eurasian-North  American  angular  velocity  vector  based  on  our  preferred  plate 
boundary  description  and  new  earthquake  fault  plane  solutions. 
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In  addition  to  the  Eurasian,  North  American,  and  Pacific  plates, 
several  smaller  plates  have  from  time  to  time  been  proposed. 

The  possible  existence  of  a China  plate  has  been  discussed  by  several 
authors  (Morgan,  1968,  1972;  Molnar  et  al.,  1973;  Das  and  Filson,  1974).  The 
primary  evidence  used  to  support  a separate  China  plate  is  the  Baikal  rift. 

Das  and  Filson  (1974)  postulate  a (west)  China  plate  rotating  clockwise  with 
respect  to  Eurasia  about  a pole  near  the  southern  tip  of  Lake  Baikal.  This 
would  account  for  active  extension  in  the  Baikal  rift,  and  some  of  the  other 
earthquake  source  mechanisms  in  Asia.  Morgan  (1972)  gives  a counterclockwise 
China-Eurasian  rotation  rate  of  2.4  X 10  deg/year  about  50°N,  127°E. 

Molnar  et  al.  (1973)  discount  the  utility  of  the  China  plate  concept  for 
describing  Asian  tectonics.  The  rift  zone  has  "spread"  no  more  than  a few 
tens  of  kilometers  since  the  beginning  of  the  Pliocene  (Horensov,  1969) . This 
corresponds  to  an  average  half-spreading  rate  of  one  to  two  orders  of  magnitude 
smaller  than  extension  at  a typical  mid-ocean  ridge.  Molnar  and  Tapponler  (1975) 
speculate  that  the  Baikal  rift  may  be  related  to  the  collision  of  India 
and  Eurasia. 

Two  additional  plates  have  been  postulated  in  the  region.  Minster  et  al. 
(1974)  proposed  a Bering  plate,  comprising  western  Alaska,  the  Bering  sea  and 
northeast  Asia,  to  explain  a systematic  misfit  of  slip  vectors  from  Aleutian 
and  7jril  trench  earthquakes  to  the  Pacific-North  American  rotation  pole. 

Though  intraplate  deformation  in  Alaska  is  no  doubt  occuring,  the  misfit  of 
slip  vectors  can  also  be  explained  by  propagation  effects  due  to  seismic  velocity 
heterogenetities  associated  with  subduction  of  the  Pacific  plate  (E.R.  Engdahl , 


Personal  communication,  1974) . We  do  not  consider  a Bering  plate  further  below. 
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Den  and  Hotta  (1973)  proposed  the  existence  of  an  Okhotsk  plate  during  the 
Mesozoic  and  early  Cenozoic  on  the  basis  of  structural  trends  and  orogenic 
belts  in  and  around  the  Sea  of  Okhotsk,  though  their  discussion  does  not 

require  a distinct  Okhotsk  plate  at  present. 

Several  alternative  plate  descriptions  for  northeast  Asia  are  considered 

in  Figure  4.  The  alternatives  shown  do  not  exhaust  all  possibilities  but 
include  most  of  those  commonly  proposed  or  assumed.  A shared  feature  of  all 
plate  models  is  that  the  boundary  between  the  Eurasian  and  North  American 
plates  continues  from  the  Nansen  ridge  onto  the  Siberian  continental  shelf 
in  the  Laptev  sea  and  along  the  active  seismic  belt  in  Yakutia  (Figure  2) . 
nemenltskava  and  Karasik  (1969),  Grachev  et.aj^  (1970),  and  Churkin  (1972) 
have  similarly  drawn  the  plate  boundary  though  this  region,  on  the  basis  of 
seismicity,  recent  faulting  and  Quaternary  volcanic  activity.  The  descriptions 
of  Figure  4 differ  in  how  the  Eurasian-North  American  plate  boundary  is  continued 

to  a triple  junction  with  the  Pacific  or  with  another  plate. 

Description  A in  Figure  4 includes  three  plates  (Eurasian,  North  American 
and  Pacific)  and  a Eurasian-North  American  boundary  through  Kamchatka.  Given 
the  EUA-NA  pole  of  rotation  (Table  2),  such  a hypothetical  boundary  would  be 
a convergence  zone.  The  boundary  should  be  a locus  of  thrust  faulting. 

However,  the  only  faults  in  Kamchatka  which  trend  in  a direction  approximately 
parallel  to  such  a proposed  boundary,  those  in  the  Kronoki-Krutogorova  fault 
zone,  are  strike-slip  in  character.  Description  A is  thus  unlikely. 

Description  B includes  four  plates:  Eurasian,  North  American,  China  and 
Pacific.  China  and  Eurasia  would  have  a boundary  striking  north-south  through 
Sakhalin.  Adopting  Morgan’s  (1972)  pole  of  rotation  for  CHI-EUA  gives  almost 
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pure  (right-lateral)  strike  slip  movement  on  meridional  faults  through 
Sakhalin;  a EUA-CHI  pole  near  Lake  Baikal  also  predicts  (left-lateral)  strike 
slip  motion  on  such  faults.  Neither  tectonic  evidence  nor  earthquake  source 
mechanisms  (below)  bear  this  out.  Description  B also  encounters  the  same 
difficulty  in  Kamchatka  as  does  case  A.  Description  B is  unlikely. 

Description  C includes  three  plates  (Eurasian,  North  American,  and  Pacific) 
with  the  EUA-NA  border  trending  north-south  through  Sakhalin  and  Hokkaido  to  a 
triple  junction  in  the  Japan  trench.  From  the  EUA-NA  pole  of  rotation,  thrust 
faults  (striking  roughly  north-south)  with  a small  amount  of  superposed  right- 
lateral  motion  would  be  expected  on  Sakhalin  and  Hokkaido,  in  good  agreement 
with  the  tectonic  evidence  discussed  above  and  the  earthquake  mechanisms  dis 
cussed  below.  We  thus  consider  this  description  to  be  an  acceptable  plate 
tectonic  interpretation  of  much  of  the  seismic  and  geologic  data,  through  it 
does  not  account  for  intraplate  deformation  in  northeast  U.S.S.R.,  China  or 
Alaska. 

Description  D has  four  plates  (Eurasian,  North  American,  China,  and 
Pacific)  with  a CHI-NA  boundary  through  Sakhalin  (cf . Mpr&flDj  1968,  1972). 
Morgan's  (1972)  CHI-NA  pole  of  rotation  is  at  54°N,  130°E  and  predicts 
(right-lateral)  strike  slip  motion  on  meridional  faults  in  Sakhalin,  in 
disagreement  with  observations.  While  one  may  argue  that  the  motion  of  the 
China  plate  has  large  uncertainties,  this  exercise  and  the  syntheses  of  Molnar 
et  al.  (1973)  and  Molnar  and  Tapponier  (1975)  imply  the  concept  of  a rigid 
China  plate  is  of  dubious  value , 

Description  E includes  the  existence  of  a separate  Okhotsk  plate  (after 
Den  and  Hotta,  1973),  along  with  the  Eurasian,  North  American,  and  Pacific 


plates.  This  interpretation  of  the  seismicity  and  tectonic  data  cannot  be 
excluded;  indeed  a cursory  examination  of  the  seismicity  in  Figure  2 
lends  support  to  the  notion  of  a separate  Okhotsk  plate,  though  the  limited 
seismicity  in  the  Sea  of  Okhotsk  precludes  a definitive  proof.  Also  because 
of  the  small  rates  and  sparse  data,  no  pole  of  rotation  for  OKH  relative  to  any 
other  plate  can  be  computed  so  the  hypothesis  of  an  Okhotsk  plate  has  no  predictive 
value.  Therefore  we  see  no  positive  reason  for  the  addition  of  another 
plate  when  the  data  can  be  explained  with  only  a three  plate  configuration. 

Of  all  the  possibilities  in  Figure  *,  then,  configuration  C is  preferable.  This 
description  is  one  of  the  simplest,  involving  only  three  distinct  plates.  The  des- 
criptlPn  provides  adequate  explanation  of  the  seismicity  end  tectonics  of  Sakhalin 
and  Hokkaido.  And  as  we  show  In  a later  section,  this  description  has  predlc  live 

4 

\ value  for  describing  the  slip  vectors  of  earthquakes  on  Sakhalin  and  Hokkaido. 

EARTHQUAKE  SOURCE  MECHANISMS 

As  an  aid  in  confirming  the  plate  description  discussed  in  the  preceding 
section  and  in  elucidating  the  tectonics  of  an  intracontinental  plate  boundary, 
we  have  determined  the  earthquake  source  mechanisms  for  all  earthquakes  in 
northeast  Asia  of  a sufficiently  large  magnitude  for  global  coverage.  We 
present  solutions  for  events  in  the  continental  shelf  in  the  Laptev  sea, 

j | 

Yakutia,  Sakhalin,  the  Tartar  strait,  and  Hokkaido. 

These  source  mechanisms  were  determined  by  utilizing  P-wave  first  motions, 

S-wave  polarizations,  and  Rayleigh  wave  amplitudes  (Figure  5,  Table  3).  The 
technique  for  determining  source  geometry  from  the  amplitude  of  the  vertical 
component  of  Rayleigh  waves  is  described  by  Forsyth  < 973) . We  utilized 
Rayleigh  wave  amplitudes  at  a period  of  67  seconds,  and  corrected  for 
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attenuation  with  a value  of  Q equal  to  125.  The  solution  which  was  compatible 
with  the  P-wave  data  and  had  a minimum  least  square  error  in  amplitude  was 
adopted  as  the  correct  mechanism. 

Event  1 occured  on  the  continental  shelf  in  the  Laptev  sea.  southeast 
of  the  Nansen  ridge.  P-wave  data  determined  one  nodal  plane,  and  surface  wave 
amplitudes  were  used  to  define  the  other  nodal  plane.  This  event  is  clearly 
a normal  faulting  event,  and  indicates  that  sea-floor  extension  occurs  on 
the  shelf  as  a continuation  of  the  Nansen  ridge.  This  earthquake  was  also 
studied  by  Conant  (1972),  who  obtained  non^rthogonal  nodal  planes  from  P-wave 
first  motions  (see  also  Figure  5),  an  artifact  of  the  heterogeneous  seismic 
velocity  structure  beneath  the  spreading  center  (Solomon  and  Julian,  1974). 

Svkea  (1967)  obtained  a similar  solution  from  an  earlier,  nearby  event. 

Events  2 and  3 occured  in  Yakutia.  Pi  Ison  ami  Frasier  (1972)  obtained 
a similar  solution  for  event  2.  The  location  of  aftershocks  (Belyy  et  .1., 

1971)  from  event  2 indicate,  that  the  fault  plane  for  that  earthquake 
strikes  northwest.  The  Inferred  solution  is  almost  pure  left-lateral 
strike  slip  motion.  Event  3,  studied  also  by  gtlstagllo  (unpublished  manuscript. 
Senior  essay,  Yale  University,  New  Haven,  Connecticut,  1974),  is  similar 
in  mechanism  to  event  2 but  has  slightly  different  nodal-plane  strikes  and 
is  less  -ell  constrained  by  the  data.  Because  of  the  similarity  in  these 
solutions  and  their  proximity,  it  might  be  assumed  that  the  northwest  striking 
plane  is  the  fault  plane  for  event  3 also.  7,.t.r«™  and  MlsHar.lna  (1965) 
also  list  strike-slip  fault  plane  solutions  for  two  Yakutia  earthquakes  in 
this  seismic  belt  (72"N,  127'E  and  66'N,  137”E).  but  we  are  unable  to  asses. 

the  quality  of  these  solutions. 

Event  4,  located  in  central  Sakhalin,  appears  to  be  almost  pure  thrust. 
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but  there  are  not  enough  data  to  tightly  constrain  the  solution.  Events  5-9 
occured  in  the  Tartar  strait  off  the  southwest  coast  of  Sakhalin.  Events  6-9 
are  part  of  an  aftershock  sequence  following  event  5.  For  this  reason,  we 
required  these  mechanism  solutions  to  be  generally  similar  while  still  satisfying 
the  data.  We  thus  obtain  five  mechanisms  which  are  all  predominantly  thrust 
events.  The  fault  plane  was  chosen  in  order  to  satisfy  several  criteria: 
the  plane  exhibit  some  right-lateral  motion  to  agree  with  the  field  geologic 
data  on  the  faults  in  Sakhalin,  the  plane  should  have  a strike  similar  to 
the  local  faults,  and  the  fault  should  be  In  agreement  with  the  shape  of  the 
isoseismal  contours  (Sslov^ev_et^li,  1973).  All  five  earthquakes  occured 
at  about  20  km  depth.  Thus  their  fault  plane  solutions  should  probably  agree 
with  the  strike  and  slip  of  the  surface  faults,  but  not  necessarily  the  dip, 
since  thrust  faults  are  commonly  shallower  in  dip  at  depth.  This  sequence 
of  earthquakes  has  the  first  two  mechanisms  similar  (5  and  6)  and  the  last 
three  (7  to  9)  all  similar  but  with  a slightly  different  fault  plane  from 
the  first  group.  All  had  a nearly  Identical  auxiliary  plane,  however.  It 
is  of  interest  that  McKenrle  (1970)  found  an  aftershock  sequence  in  the 
Mediterranean  region  which  similarly  had  a constant  slip  vector,  but  a 
differing  fault  plane,  for  each  individual  event. 

Event  10  occured  In  Hokkaido  at  a depth  of  25  km  beneath  the  Hidaka 
mountains.  We  consider  it  to  be  unrelated  to  underthrusting  of  the  Pacific 
plate  because  of  its  shallow  depth  and  even  shallower  aftershock  sequence,  and 
because  of  the  orientation  of  the  fault  plane.  From  the  aftershock  distribution 
fMorlva,  1972)  the  shallowly  dipping  plane  was  determined  to  be  the  fault 
plane.  It  thus  was  a thrust-fault  event  with  a component  of  left-lateral 

Btrike  slip  motion. 


A logical  question  is  whether  these  earthquake  mechanism  solutions  are 
compatible  with  the  earlier  discussion  of  plate  boundaries  and  with  the  EUA-NA  j 

pole  of  rotation.  Event  1 implies  that  such  a vole  must  be  located  south  of 
76.5°N,  in  order  to  have  extension  on  that  region  of  the  shelf.  If  events  2 


! i 


and  3 are  both  Interpreted  as  left-lateral  faulting  occuring  on  a single 
plate  boundary,  the  two  slip  vectors  uniquely  define  a rotation  pole  at  65'N, 
148”E.  Such  a pole,  however,  would  be  in  systematic  disagreement  with  fracture 
rone  trends  and  earthquake  slip  vector  data  in  the  Atlantic  and  Arctic  oceans. 
The  pole  would  be  well  outside  the  951  confidence  ellipse  of  Minster  et  .L 
(1974)  . Consequently  these  two  events  cannot  be  indicative  of  rigid  motion  of 
the  Eurasian  and  North  American  plates.  Rather  they  are  the  product  of 
complicated  deformation  in  a boundary  between  two  plates  near  the  relative 
rotation  pole  for  the  same  plates.  He  comment  further  on  this  point  in  the 
following  section. 

If  the  plate  tectonic  description  of  northeast  Asia  is  as  discussed 
above  (case  C,  Figure  5),  then  sufficiently  far  from  the  EUA-NA  pole,  the  slip 
vectors  of  earthquakes  on  the  boundary  should  be  predictable.  Using  the 
EUA-NA  angular  velocity  vector  of  Minster  et  al^  (1974),  the  slip  vectors 
for  earthquakes  on  Sakhalin  and  Hokkaido  should  have  azimuths  of  about  75“ 
to  80°.  These  values  are  very  close  to  those  observed  (Figure  5,  Table  3), 
lending  substantial  credence  to  the  above  identification  of  the  North  American- 
Eurasian  plate  boundary. 

A logical  next  step  is  to  recalculate  the  EUA-NA  rotation  pole  using  the 
new  data.  Combining  the  slip  vectors  for  events  4,  5 and  7 through  10  (the 
slip  vector  for  event  6 is  poorly  constrained)  with  essentially  the  data  set 
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(Table  4)  for  EUA-NA  rotation  of  Minster  et  "l.  (1974),  we  obtain  a EUA-NA 
pole  at  61.8°N,  130. 0°E,  with  a rate  of  2.48  X 10~?  deg/yi  (Figure  7).  Such 
a location  is  approximately  7°  south  of  the  Minster  et  al. (1974)  pole  but  is 
within  their  95%  confidence  ellipse.  Therefore  it  is  not  a statistically 
significant  improvement  over  the  Minster  et  al.  (1974)  pole,  rather  it  is  a 
pole  of  rotation  which  explains  a larger  data  set.  This  pole  of  rotation  describes 
the  relative  motion  of  the  Eurasian  and  North  American  plates  in  the  Atlantic, 
the  Arctic,  and  in  Sakhalin  and  Hokkaido.  It  does  not  describe  the  moticn 
within  about  10°  of  the  rotation  pole  (e.g.,  events  2 and  3). 

ON  THE  NATURE  OF  INTRACONTINENTAL  PLATE  BOUNDARIES 

The  intracontinental  portion  of  the  boundary  between  the  Eurasian 
and  North  American  plates  has  several  characteristics  which  distinguish  it 
from  the  more  common  submarine  plate  boundaries.  We  comment  on  these 
characteristic;1  in  this  section,  with  occasional  generalizations  to  other 
intracontinental  plate  edges. 

The  plate  boundary  in  Yakutia  is  very  wide  and  diffuse.  At  its 
widest,  the  boundary  (if  indeed  such  a term  is  still  appropriate)  is  600 
km  wide.  The  diffuse  nature  of  the  boundary  is  more  likely  a property  of 
continental  lithosphere  than  duo  to  the  slow  relative  plate  velocity.  The 
width  of  the  seismic  zone  increases  markedly  between  the  Nansen  ridge  (oceanic 
lithosphere)  and  its  extension  onto  the  continental  shelf  in  the  Laptev 
sea  (Figures  1 and  2) . Other  intracontinental  plate  boundaries  of  different 
types  and  with  different  relative  plate  velocities  share  this  very  extended 
character;  western  North  America  (Atwater,  1970)  is  a good  example.  That 
continental  lithosphere  is  generally  thicker  than  oceanic  lithosphere  may 
be  part  of  the  answer  for  the  diffuse  definition  of  the  intracontinental 
edges  of  plates.  More  important,  probably,  is  that  continental  lithosphere 
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is  very  heterogeneous,  a complex  cementation  of  blocks  and  belts  of  different 
make-up,  texture  and  age.  Compared  to  the  relatively  fresh  and  relatively 
homogeneous  oceanic  plates,  continents  have  generally  undergone  a long  history 
of  stress-  and  fracture-producing  tectonic  activity  and  are  crisscrossed  with 
weak  zones  highly  susceptible  to  deformation  when  stressed. 

The  relative  rotation  pole  is  near  the  plate  boundary . It  is  sometimes 
difficult  to  separate  cause  and  effect  in  a physical  phenomena,  but  we 
speculate  that  the  location  of  the  Eurasian-North  American  rotation  pole 
very  near  the  intracontinental  boundary  between  these  two  plates  is  not 
coincidence.  Rather  it  is  likely  related  to  a greater  asthenospheric 
resistance  to  moving  continental  than  oceanic  lithosphere  in  general 
(Solomon  et  al.,  1975;  Forsyth  and  Uyeda,  1975)  and  to  a difficulty  in 


subducting  one  continental  block  beneath  another  (McKenzie,  1969),  particularly 
without  first  closing  an  intervening  ocean.  The  situation  is  not  quite  that 
simple,  since  the  plates  are  responding  to  a number  of  different  types  of 
forces  (Solomon  and  Sleep,  1974) . Nonetheless  there  are  at  least  two  other 


possible  examples  cf  where  a relative  rotation  pole  is  located  near  the 
intracontinental  portion  of  a plate  boundary:  the  Pacific-Indian  pole  is  not 
far  from  New  Zealand  and  the  African-Eurasian  pole  is  not  far  from  the 
straits  of  Gibraltar. 

The  displacements  during  large  earthquakes  are  not  predictable  near  the 
relative  rotation  pole.  This  characteristic  is  related  to  the  p evious  two. 
Because  the  rotation  pole  is  near  the  plate  boundary,  the  stress  system  in  the 
boundary  zone  changes  rapidly  with  small  changes  in  distance.  Structural 
heterogeneities  modulate  the  stress  field  and  the  materiel  response  to  stress. 
Aggravating  the  complexity  of  the  stress  and  stress  release  fields  is  that  the 
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ins  tan  taneous  rotation  pole  is  not  fixed,  but  migrates  with  respect  to  the 
two  plates  (Pitman  and  Taiwan!.  1972) . 

The  modern  boundary  is  closely  related  to  ancient  plate  boundaries.  The 
broad  seismic  belt  in  Yakutia  (Figure  2)  marking  the  location  of  the  current 
Eurasian-North  American  plate  boundary  lies  within  the  Cherskiy-Verkhoyansk 
foldbelts,  which  Churkin  (1972)  has  interpreted  as  a fossil  suture  marking 
the  early  Cretaceous  collision  of  two  continental  blocks • The  current  plate 
boundary  through  Sakhalin  and  Hokkaido  (Figure  6)  follows  closely  a Mesozoic 
plate  boundary  marking  the  locus  of  eastward  subduction  of  one  plate  beneath 
another  (Sugimura  and  Uyeda,  1973;  Den  and  Hotta,  L973) . This  is  a familiar 
story  in  the  plate  tectonic  evolution  of  the  earth's  surface:  fossil  plate 

boundaries  are  apparently  relatively  weak  portions  of  continental  blocks  and 
are  the  preferred  sites  for  cveation  of  new  plate  edges. 


CONCLUSIONS 

Though  plate  boundaries  within  continents  can  rarely  be  defined  with 
precision,  the  boundary  between  the  Eurasian  and  North  American  plates  between 
the  Arctic  and  Pacific  oceans  can  be  identified  on  the  basis  of  seismicity, 
recent  tectonics  and  earthquake  source  mechanif ms . The  simplest  plate 
configuration  that  adequately  accounts  for  these  data  continues  the  Eurasian- 
North  American  boundary  from  the  Nansen  ridge  through  a broad  seismically 
active  zo:>e  in  northeast  U.S.S.R.  to  the  Sea  of  Okhotsk  and  thence  southward 
through  Sakhalin  and  Hokkaido.  With  this  configuration,  the  slip  vectors 
derived  from  earthquake  mechanisms  in  Sakhalin  and  Hokkaido  are  predicted  from 
the  EUA-NA  relative  rotation  pole.  A new  pole,  only  slightly  different  from 
other  recent  solutions,  is  computed  on  the  basis  of  the  additional  sl'p  vector 
data. 


i£ *... 


-80- 


I 


1 


The  intracontinental  plate  boundary  is  spread  over  a width  of  as  much  as 
600  km,  and  tends  to  follow  ancient  plate  margins.  Deformation  in  the  vicinity 
of  the  rotation  pole,  which  lies  near  the  boundary,  is  poorly  described  by 
rigid-plate  tectonics.  These  features  can  be  explained  by  the  thickness  anc 
heterogeneity  of  continental  lithosphere  and  by  the  influence  of  continents  on 
the  forces  moving  the  plates. 
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TABLE  1 

RELOCATED  EARTHQUAKES  IN  NORTHEAST  USSR 


Date 

Origin  Time 
(G.M.T.) 

Latitude 

Longitude 

h 

m 

s 

Nov.  30,  1918 

06 

48 

38.2 

70.704°N 

133 . 363°E 

Mar.  13,  1924 

10 

41 

58.7 

62.772 

150.062 

Mar.  15,  1924 

10 

31 

21.3 

49.176 

142.570 

May  27,  1924 

20 

09 

30.3 

62.452 

135.056 

Feb.  18,  1925 

11 

36 

3.7 

66.614 

145.648 

Apr.  9,  1926 

10 

4 

32.0 

72.865 

132.093 

June  10,  1927 

18 

13 

23.4 

48.364 

139.067 

Nov.  14,  1927 

00 

12 

7.4 

70.233 

128.733 

Nov.  14,  1927 

04 

56 

29.5 

70.208 

128.990 

Nov.  15,  1927 

21 

48 

45.7 

70.275 

129.069 

Feb.  3,  1928 

13 

47 

36.5 

70.374 

128.126 

Feb.  21,  1928 

19 

49 

6.0 

67.573 

-172.561 

Feb.  24,  1928 

14 

10 

25.4 

67.536 

-173.824 

Feb.  26,  1928 

01 

19 

12.8 

67.395 

-171.034 

Aug.  16,  1928 

07 

36 

44.6 

69.842 

123.130 

Aug.  25,  1928 

01 

48 

32.1 

49.060 

141.814 

July  15,  1931 

16 

27 

0.6 

59.082 

148.185 

Oct.  10,  1931 

16 

37 

8.4 

59.504 

148.027 

July  10,  1932 

00 

43 

26.3 

52.642 

142.052 

Sept.  7,  1933 

22 

39 

20.3 

61.963 

177.429 

Oct.  25,  1935 

17 

38 

14.6 

51.854 

142.887 

Nov.  3,  1936 

04 

43 

23.2 

59.198 

152.815 

May  10,  1947 

00 

07 

14.5 

57.858 

141.908 

Apr.  14,  1951 

13 

33 

2.1 

61.117 

136.306 

Note:  All  depths  were  constrained  to  be  15  km 
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FIGURE  CAPTIONS 

Figure  1.  Seismicity  of  the  Arctic,  1962-1969,  after  Nat.  Ocean.  Atmos. 

Adm.  (1970) . 

Figure  2.  Shallow  seismicity  of  northeast  Asia,  1909-1973,  azimuthal  equal-area 

projection  about  69.3°N,  128°E.  The  size  of  a symbol  is  proportional  to 
(body-wave)  magnitude  M.  Triangles  indicate  a seismograph  station. 

Only  events  between  latitudes  41°  and  75°N  and  between  lotgitudes 
120°  and  170°E  are  included.  For  clarity,  no  events  are  plotted  in  the 
Aleutian  and  Kuril  trenches  with  M < 5.5  or  on  Sakhalin  with  M < 4.0. 

Figure  3.  Cenozoic  tectonics  of  the  Sea  of  Okhotsk  region,  simplified  from 
Alverson  et  al.  (1967) . 

Figure  4.  Possible  plate  descriptions  of  the  seismicity  and  tectonics  in 

northeast  Asia.  The  dashed  lines  indicate  boundaries  between  plates. 
EUA  - Eurasian,  NA  - North  American,  PAC  - Pacific,  CHI  - China,  OKH  - 
Okhotsk. 

Figure  5.  Earthquake  source  mechanisms.  Numbers  refer  to  the  event  number  in 
Table  2.  Equal  area  projection  of  P-vave  first  motions  and  S-wave 
polarizations  plotted  on  the  lower  focal  hemisphere;  closed  circles 
are  compression,  open  circles  are  dilatation.  Triangles  indicate 
short  period  data  from  the  Bulletin  of  the  Japanese  Heterological 
Agency  (event  4 only) . All  other  data  were  read  from  long  period 
WWSSN,  Canadian  network,  and  Lamont  network  seismograms.  Rayleigh- 
wave  amplitude  radiation  patterns  are  shown  for  events  1 and  6. 

Figure  6.  The  Eurasian-North  American  plate  boundary  in  northeast  Asia.  The 
newly  computed  EUA-NA  pole  of  rotation  is  shown  together  with  the 
95Z  confidence  ellipse.  Locations  of  the  pole  of  rotation  from 
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4.  THREE  DIMENSIONAL  EARTH  STRUCTURE 

t 

4.1  Determination  of  the  Three  Dimensional  Structure  of  the 
Lithosphere  by  K.  Aki,  A.  Christof fersson  and  E.  Husebye 
A new  three-dimensional  earth  modelling  is  formulated 
in  order  to  meet  an  increasing  demand  for  more  detailed  and 
accurate  information  about  the  earth's  interior.  We  start 
with  a layered  medium  of  classic  seismology,  but  divide  each 

I 

layer  into  many  blocks  and  assign  a parameter  to  each  block 
which  describes  the  velocity  perturbation  from  the  average 
for  the  layer.  Our  data  is  the  teleseismic  P travel  time 
t residuals  observed  at  an  array  of  seismographs  distributed 

on  the  surface  above  the  earth's  volume  we  are  modelling. 

By  isolating  various  sources  of  errors  and  biases,  we  arrive 
at  a system  of  equations  to  determine  the  model  parameters. 

The  solution  was  obtained  by  the  use  of  generalized  inverse 
and  stochastic  inverse  methods  with  the  analysis  of  resolution 
and  errors  in  estimates.  Our  method  also  gives  a lower  limit 
of  the  true  r.m.s.  velocity  fluctuation  in  the  actual  earth 
under  the  assumption  of  ray-theory. 

Using  NORSAR  P-wave  residuals,  we  have  determined  the 
three-dimensional  seismic  structure  of  the  lithosphere  under 
the  array  to  the  depth  of  126  km.  The  true  r.m.s.  velocity 
fluctuation  was  found  to  be  at  least  3.4%.  This  is  in  agree- 
ment with  estimates  obtained  from  statistical  analysis  of 


r 


* 


P time  fluctuation  based  on  the  Chernov  theory.  The  three- 
* dimensional  velocity  anomalies  are  presented  both  by  the 

generalized  inverse  and  by  the  stochastic  inverse  solutions. 
We  preferred  the  dual  presentation,  because  it  gives  the 
reader  greater  freedom  in  judging  the  results  than  a single 
"optimal"  solution.  Both  methods  gave  essentially  the  same 
results.  The  discrepancies,  when  they  existed,  were  always 
explainable  in  terms  of  differences  in  the  smoothing  proce- 
dure which  is  explicitly  given  in  the  resolution  matrix. 

We  found  clear  evidence  of  ^ ipe-like  structures  under 
NORSAR,  dipping  northward  and  away  from  the  surface  contour 
of  the  Oslo  graben.  These  pipe-like  structures  were  inter- 
preted as  vestiges  of  magma  ascent  by  penetrative  convections 
associated  with  the  Permian  volcanism  of  the  Oslo  graben. 

The  inclination  of  the  pipe-like  structures  is  interpreted 
as  a result  of  plastic  deformation  of  the  lithosphere  due 
j to  the  shear  exerted  by  the  asthenosphere  convection  current 

driving  trie  plate  motion. 
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4 . 3 Three-Dimensional  Seismic  Structure  of  the  Lithosphere 
under  Montana  LASA  by  K,  Aki,  A.  Christoff ersson  and 
E.  Husebye 
ABSTRACT 

Using  P-wave  residuals  for  teleseismic  events  observed 
at  the  Montana  LASA,  we  have  determined  the  three- dimensional 
seismic  structure  of  the  lithosphere  under  the  array  to  a 
depth  of  140  km.  The  r.m.s.  velocity  fluctuation  was  found 
to  be  at  least  3.2%  which  may  be  compared  to  estimates  of  ca 
2%  based  on  the  Chernov  random  medium  theory.  The  solutions 
are  given  by  both  the  generalized  inverse  and  stochastic 
inverse  methods  in  order  to  demonstrate  the  relative  merit 
of  different  inversion  techniques. 

The  most  conspicuous  feature  of  the  lithosphere  under 
LASA  is  a low  velocity  anomaly  in  the  central  and  NE  part 
of  the  array  siting  area  with  the  N60°E  trend  and  persisting 
from  the  upper  crust  to  depths  greater  than  100  km.  We 
interpret  this  low  velocity  anomaly  as  a zone  of  weakness 
caused  by  faulting  and  shearing  associated  with  the  building 
of  the  Rocky  Mountains. 
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Introduction 


The  seismic  velocity  structure  beneath  LASA  has  intrigued 
seismologists  ever  since  this  large  aperture  array  was 
established  in  Montana  in  1965  The  interpretation  of  two 
dimensional  seismological  observations  was  novel  to  seismology 
and  resulted  in  numerous  papers  dealing  with  the  crust  and 


upper  mantle  under  LASA  For  references , see  Toksftz,  Chinnery 
and  Anderson  (1067);  Mack  (1969)/  Greenfield  and  Sheppard 
(1969)/  Glover  and  Alexander  (1969)/  Larner  (1970)/  Engdanl 
and  Felix  (1971)/  Iyer  and  Healy  (1972);  Warren,  Healy,  Bohn 
and  Marshall  (1973);  Aki  (1973);  and  Capon  (1974)  Besides 
partly  conflicting  results  these  papers  are  facinating 
reading  as  they  reflect  an  area  of  progress  in  seismology 


but  at  the  same  time  serve  to  expose  the  enormity  of  the  non- 


uniqueness problem  inherent  to  inverting  seismic  observations. 


Experience  with  LASA  has  shown  that  conventional  layered 
medium  models  used  for  interpreting  seismic  refraction  and 
wide-angle  reflect: c1'  data  cannot  satisfactorily  explain  the 
observed  P-.wave  travel  time  and  amplitude  anomalies. 

Two  fundamental  problems  are  encountered  in  inversion  of 
teleseismic  P-wave  travel  time  data*  namely  spatial  resolution 
and  flexibility  in  the  basic  model  formulation.  In  the  former 


case,  we  are  concerned  with  the  location  of  the  structural 


i 


l 


IS 


V 


’ 


I 


inhomogeneities,  e,g,,  whether  the  observed  time  residuals  principally 


I 
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originate  in  the  lithosphere  beneath  the  array  (Iyer  and 
Healy,  1972),  in  the  mantle  (Toksfiz  et  al . f 1967',  in  the 
lowermost  part  of  the  mantle  or  in  the  source  region  (Davies 
and  Sheppard,  1972).  In  the  case  of  model  formulation,  it 

I 

should  be  pointed  out  that  conventional  layered  models  with 
uniform  material  properties  cannot  accomodate  complex 
geological  features.  In  fact,  the  antipode  to  conventional 
approaches  is  to  assume  that  the  seismic  structure,  say  of  the 
lithosphere,  is  so  complex  that  it  behaves  as  a medium  with 
random  inhomogeneities  This  type  of  model  has  successfully 
been  introduced  in  the  interpretation  of  LASA  P'WaVw  travel 
time  and  amplitude  anomalies,  by  using  the  Chernov  (1960) 
theory  for  wave  propagation  in  random  media  (eg.,  see  Aki 
197 3 y Capon,  1974;  and  Berteussen,  Christoffersson,  Husebye 
and  Dahle,  1975) . The  Chernov  theory  is  based  on  a statistical 
description  of  the  medium  in  terms  of  the  mean  square  of  fractional 

fluctuation  of  velocity,  the  correlation  distance  (the  distance 

i 

between  two  points  in  the  medium  at  which  velocity  fluctuation 
becomes  approximately  uncorrelated) f and  the  linear  extent  of 
the  inhomogeneous  region. 

In  this  paper,  we  shall  utilize  our  new  earth  modelling 
approach  with  velocity  fluctuations  similar  to  the  Chernov 
medium  but  defined  deterministically  within  the  volume  of 

.• 
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crust  and  upper  mantle  for  estimating  the  three-tdimensi cnal 
seismic  structures  beneath  LASA  First,  we  shall  give  a 
brief  description  of  our  earth  modelling  technique  and  then 
discuss  our  results  in  relation  to  available  geophysical 
and  geological  information  pertinent  to  the  area, 

Model  Formulation 

In  this  section  a brief  presentation  is  given  of  our 
earth  modelling  approach  and  the  inversion  technique  of 
P-wave  travel  time  residuals  used  for  estimating  velocity 
perturbations  in  crust  and  upper  mantle  (see  Aki,  Christoff erton 
and  Husebye,  1975,  hereafter  referred  to  as  Paper  1).  The 
Initial  model,  shown  in  Fig.  1,  consists  of  homogeneous  layers 
with  constant  thickness  and  fixed,  average  velocity  and 
comprises  a somewhat  arbitrarily  defined  volume  under  the  array 
containing  the  observed  ray  paths  The  earth’s  structure  outside 
the  model  (see  Fig,  1}  is  denoted  ‘standard  earth ,f  and  is 
assumed  known.  Each  layer  of  the  initial  model  is  divided  into 
many  blocks,  and  within  each  layer  the  relative  velocity 
perturbation  m for  a given  block  are  defined  as  the  average 
perturbation  over  all  rays  having  most  of  the  unperturbed  path 
within  the  block. 

Omitting  all  details  which  are  given  in  Paper  1,  the  model 
for  estimating  velocity  perturbation  for  the  individual 
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blocks  is  given  by: 


T = G m + t 


(1) 


where  the  elements  of  the  T -vector  are  the  travel  time 
residuals  minus  their  average  over  all  sensors  for  each  event. 
The  G-matrix  are  calculated  from  the  travel  time  each  ray 
spends  within  each  layer  for  the  initial  unperturbed  model. 

The  m-vector  contains  the  unknown  block  velocity  perturbations 
and  the  e-vector  the  higher  order  and  error  terms.  The  sum 
of  the  columns  in  G corresponding  to  a layer  is  identically 
zero,  and  therefore  adding  a constant  to  the  velocity 
perturbations  within  a layer  cannot  affect  T.  Thus  the  rank 
of  G is  less  or  equal  to  the  number  of  blocks  sampled  by  the 
observational  data  (or  hit  by  the  incoming  rays)  minus  the 
number  of  layers  in  the  model.  This  means  that  the  velocity 


i 


Using  the  method  of  generalized  inverse  (Lancsoz,  1961), 
the  solution  is: 


, T _g  t 

« (G  G)  * G T (3) 

The  resolution  matrix  R (Bachus  and  Gilbert,  1968)  is  given 
by 

R ■ WT  (4) 

where  the  V-matrix  contains  the  eigenvectors  corresponding 
to  the  non-zero  eigenvalues  of  the  (GTG) -matrix.  In  case 
of  the  generalized  inverse  solution  the  elements  of  R 
corresponding  to  the  best  possible  resolution,  are: 
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in 

different  layers 

(n  - number  of  blocks  sampled  in  a layer) 

If  we  assume  that  the  covariance  of  the  residual  time  or 
the  £- vector  in  Equation  (1)  can  be  written  as  o2I,  the 
standard  errors  of  the  mg-estimate  would  be: 

[D2  (mg)] 1/2  = (oV'V) 1/2  (5) 

where  the  A -matrix  contains  the  non-zero  eigenvalues  of 

r 

T 

the  (G  G) -matrix.  It  is  shown  in  Paper  1 that 
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(U  ) 
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NB 


mg 


2 1/2 
- tr  D (mg ) ) ] 


(6) 


represents  an  estimate  of  a lower  limit  for  the  true  root  mean 

1/2 

square  slowness  perturbation  (y^)  in  the  volume  of  the 
earth  sampled.  NB  is  the  total  number  of  blocks  sampled 
and  tr  represents  the  trace  of  a matrix. 

An  alternative  to  the  generalized  inverse  approach  is 
a special  case  of  the  stochastic  inverse  solution  (Franklin, 
1970) , namely: 

S^st  = Wtg  + e2!) 'Vt  (7) 

This  solution  has  smaller  standard  errors  but  poorer 
resolution  than  the  generalized  inverse  as  discussed  in 
Paper  1.  In  particular,  r^j  is  no  longer  zero  for  bloc} s in 
different  layers.  This  will  introduce  vertical  smoothing 
between  blocks  in  different  layers  which  in  turn 
introduce  distortion  in  the  solution  as  the  smoothing  kernels 
varies  from  block  to  block.  Both  solutions  have  been  calculated 
for  LASA  and  will  be  presented  in  the  next  section. 

Finally,  let  us  discuss  possible  shortcomings  and  biases 
in  our  earth-modelling  approach  and  consequently  in  the 
estimated  velocity  perturbation  solution.  As  ray-theory  is 
used  for  computing  the  travel  time  through  the  initial  model 
(see  Fig.  1),  it  is  implicit  that  the  inhomogeneities  are 
smooth  within  a wavelength.  On  the  other  hand,  possible 
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inhomogeneities  with  scale  length  much  shorter  than  the 
wavelength  will  behave  as  an  equivalent  homogeneous  body  with 
some  average  properties.  Thus,  it  is  meaningless  to  make 
the  block  size  much  smaller  than  a wavelength.  Our  block 
size  of  20  km,  however,  is  determined  primarily  by  the  station 
density.  Ideally,  we  would  like  to  divide  the  whole  earth 
into  blocks,  but  the  finite  amount  of  data  available  imposes 
limits  to  the  number  of  blocks  and  volume  of  the  earth  to  be 
included  in  our  model.  Because  of  this  limitation,  the  effect 
of  possible  deeper  inhomogeneities  may  be  projected  into  our 
model  thus  biasing  the  final  results.  Besides  using  many  events 
having  a reasonable  azimuth-distance  coverage,  two  quantitative 
measures  for  checking  such  biases  are  discussed  in  Paper  1. 
First,  extending  our  model  from  4 to  5 layers  does  not 
significantly  change  the  velocity  perturbations  in  the  4 
uppermost  layers;  secondly  the  model  explains  about  85%  of 
the  variance  in  the  observed  travel  time  residuals,  and  the 
remaining,  unexplained  variance  is  comparable  to  the  variance 


of  reading  error. 

Inversion  of  LAS A P-wave  Residuals 

The  LASA  array  in  Montana  consists  of  21  subarrays  and 
its  configuration  is  shown  in  Figure  2.  Travel  time  residuals 
of  short  period  P-waves  at  the  subarray  centers  (with  respect 
to  Herrin  tables)  for  a large  number  of  teleseismic  events 
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have  been  tabulated  by  Chiburis  and  Ahner  Cl 9 7 3 ) . This 
data  set,  used  as  station  corrections  for  beam  steering  during 
routine  array  processing,  is  based  on  carefully  analyzed  events 
and  the  measurement  errors  are  considered  to  be  less  than  Q.l  sec. 
Moreover,  the  above  events  have  an  adequate  azimuthal  and 
distance  coverage  (see  figure  3)  so  these  observations  are  well 
suited  for  our  inversion  experiment.  The  LASA  data  used 
comprises  178  events  recorded  at  17  subarrays  giving  a total 
of  3.026  readings.  The  4 F-ring  seismometers  were  excluded 
from  analysis  due  to  their  relatively  large  separation  from 
the  other  stations.  Exa.  pies  of  observed  travel  time  residuals 
are  listed  in  Table  1.  The  initial  model  parameters  for  the 
crust  and  upper  mantle  beneath  LASA  are  given  in  Table  2,  and 
is  based  on  the  results  of  Warren  et  al . , (1973)  and  Tcksttz 
et  al.  (1967).  Each  of  the  five  layers  of  the  initial  model  is 
divided  into  9x9  square  blocks  with  side  length  20  km. 

We  shall  present  both  the  generalized  c.nd  stochastic 
inverse  solutions.  This  dual  presentation  of  our  results  is 
preferred  to  a single  optimal' solution  based 
on  the  trade-off  between  resolution  and  standard  errors  as 
discussed  by  Backus  and  Gilbert  (1970)  . Moreover,  as  the 
amount  of  data  available  is  limited,  choice  of  block  configuration 
may  affect,  the  cross-sampling  of  the  individual  blocks,  and 
consequently  the  final  solution.  One  way  of  avoiding  this 
effect  is  to  average  the  results  for  say  different  block 
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configurations.  This  block  averaging  procedure  was  adopted  for 
the  stochastic  inverse  solution,  but  not  for  the  generalized 

: 

inverse  in  order  to  maintain  the  best  possible  resolution. 

The  generalized  inverse  solution  is  given  in  Figures  4-8 
and  also  tabulated  together  with  estimated  standard  errors  in 
Table  3.  The  velocity  perturbations  are  expressed  in  percentage 
of  average  layer  velocity  (given  in  Table  2) , and  relatively 
high  (H)  and  low  (L)  velocities  are  given  negative  and  positive 
signs  respectively.  Shaded  blocks  represent  significantly 
high  and  low  velocity  areas  with  the  magnitude  of  velocity 
perturbations  exceeding  twice  the  corresponding  standard  errors. 

All  but  five  of  the  sampled  blocks  are  perfectly  resolved  in 
the  generalized  inverse  solution  except  for  the  uncertainty  of 
an  additive  constant  for  each  layer  as  discussed  in  the 
preceeding  section.  For  all  blocks  except  the  5 encircled  in 
Figures  4-8,  the  elements  of  the  calculated  resolution  matrix 
are  in  accordance  with  Equation  (3) , demonstrating  that  the 
resolution  is  the  best  possible.  The  numerical  values  of  the 
diagonal  elements  were  0.96-0.99  and  the  off-diagonal  elements 
were  * (0. 01-0.04) . Elements  corresponding  to  combinations  of 
blocks  in  different  layers  were  zero.  The  part  of  the 
resolution  matrix  corresponding  to  the  5 blocks  not  perfectly 
resolved,  is  listed  in  Table  4:  the  overall  fit  of  the 
generalized  inverse  solution  to  the  observational  data  is 

I 

excellent.  The  standard  deviation  of  the  residuals  is  0.08  sec 
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which  is  comparable  to  the  measurement  error*  The  estimate 
of  the  lower  limit  of  the  root  mean  square  of  the  true 
velocity  perturbations  is  3.2%  under  LASA  as  compared  to 
3.4%  under  NORSAR, 

The  stochastic  inverse  was  calculated,  with  the  damping 
factor  0^  = 0.02  (sec  per  %1^,  for  two  sets  of  block  con- 
figurations, which  were  identical  except  for  the  shift  in  the 

2 

NE-SW  directions  by  a half  block  size.  This  value  of  6 is 
125  times  greater  than  the  smallest  non-zero  eigenvalue  of 

0. 00016.  Out  of  the  total  271  non-zero  eigenvalues,  104  are 
smaller  than  0 , This  gives  a rough  idea  about  the  resolution 
and  standard  errors  in  the  stochastic  inverse  solution  through 
Equation  (49)  and  (50)  in  Paper  1.  The  final  solution  was 
obtained  by  the  averaging  over  four  partly  overlapping  blocks, 

1. e.  two  from  the  'unshifted ' solution  and  two  from  the  'shifted' 
solution.  The  averaged  value  is  plotted  at  the  center  of  gravity 
of  the  four  blocks  and  subsequently  interpolated  anomaly  contours 
at  1%  intervals  are  drawn  as  shown  in  Figs.  9-13. 

An  example  of  the  resolution  matrix  for  the  stochastic 
inverse  is  shown  in  Fig.  14.  On  the  extreme  left,  there  are 
five  tic-tac-toe  meshes  corresponding  to  the  five  layers  of 
our  model.  The  center  of  each  mesh  corresponds  to  the  block 
at  a geographic  location  "A"  shown  in  the  map  at  the  right 
hand  of  the  same  figure.  These  five  meshes  show  elements  of  the  row 
vector  of  the  resolution  matrix,  for  which  the  diagonal  element 
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corresponds  to  the  block  in  Layer  1 at  geographic  location  ’A" . 
The  figure  shows  that  the  diagonal  element  is  84%,  and  the 
element  corresponding  to  the  blocks  in  the  same  layer  as  the 
diagonal  element  show  small  negative  values  as  also  observed 
in  the  case  of  the  generalized  inverse  solution.  The  elements 
for  blocks  immediately  below  the  diagonal  element  block 
exhibit  a small  positive  value  of  5%  indicating  a slight 
vertical  smoothing  for  block  "A"  of  Layer  1.  The  elements 
not  shown  in  the  figure  are  all  less  than  1%.  Figure  14  also 
shows  elements  of  the  resolution  matrix  for  four  other  blocks 
with  geographic  locations  indicated  by  alphabetic  letters  and 
depth  indicated  by  layer  number.  Notice  that  the  resolution 
is  very  good  for  the  central  blocks  (except  for  Layer  5)  but 
may  be  poor  for  spherical  blocks  as  demonstrated  in  Table  5. 

In  the  latter  case  vertical  smoothing  effects  are  significant. 


Interpretation  of  3-dimensional  seismic  velocity  anomalies  under  LASf 


The  seismic  network  of  LASA  is  located  in  the  transition 
zone  between  the  Rocky  Mountains  and  the  Great  Plains.  Since 
the  Rocky  Mountains  trend  roughly  N-S,  one  might  expect  a gradual 
change  of  the  crust  and  upper  mantle  in  the  E-W  direction.  Our 
results,  however,  demonstrate  that  changes  are  neither  gradual 
nor  trending  in  the  E-W  direction. 

Before  describing  the  velocity  anomalies  presented  in  the 
previous  section,  we  shall  summarize  available  geological  and 
geophysical  data  pertinent  to  the  area.  The  LASA  siting  area 
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is  covered  by  undisturbed  sediments  of  thickness  about  3 km 
lying  conformably  on  Precambrian  basement  rocks  (Brown  and 
Poort,  1965).  The  isopach  lines  for  the  depth  to  basement 
trend  roughly  NW-SE,  increasing  toward  east.  Glover  and 
Alexander  (1969)  suggested  the  same  trend  for  deeper  structure 
using  limited  data  on  Rayleigh  wave  dispersion  and  spectral 
ratio  for  long  period  P waves,  but  they  also  recognized  the 
NE-SW  trend  dominant  in  the  travel  time  and  amplitude  anomalies 
of  short  period  P waves.  The  spectral  ratio  method  was  used  in 
analysis  of  LASA  data  also  by  Bakun  (1971)  whose  results  showed 
a local  NE-SW  trend  near  the  array  center  and  a regional  NW-SE 
trend.  Engdahl  and  Felix  (1971)  convincingly  demonstrated  that 
the  main  source  of  the  travel  time  anomalies  lie  in  the  crust  and 
upper  mantle  beneath  the  array  by  showing  that  practically 
identical  anomaly  patterns  prevail  for  waves  having  greatly 
different  paths  through  the  deep  interior  of  the  earth  but 
sharing  nearly  common  ray  paths  near  the  array. 

In  an  early  attempt  to  explain  the  travel  time  residuals, 
Greenfield  and  Sheppard  (1969)  proposed  a model,  in  which  the 
crust  thickens  sharply  from  58  to  70  km  toward  NW  near  the 
center  of  array.  This  two-dimensional  model  is  trending  N60°E, 
and  the  crust  is  thinning  gradually  both  toward  N30°W  and  S30°E 
and  away  from  the  center  of  array,  A similar  model  was  postulated 
by  Larner  (1970) , who  applied  the  wave-theoretical  method  of  Aki 
and  Larner  (1970)  for  LASA  structural  analysis  based  on  the  short 
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period  P-wave  travel  time  and  amplitude  observations.  Using 
additional  travel  time  data  from  the  extended  U.S.  Geological 
Survey  array,  Iyer  and  Healy  (1972)  presented  a detailed  map 
of  the  crustal  thickness  under  LASA.  Their  result  essentially 
agrees  with  Larner's  i.e.,  the  crust  has  a maximum  thickness 
of  56  km  under  the  array  center  but  is  thinning  to  about  40  km 
both  to  the  north  and  to  the  south. 

The  E-NE  trending  crustal  structure  variation  under  LAS.i 
has  gained  support  'rom  numerous  seismic  refraction  and  wide 
angle  reflection  studies,  summarized  by  Warren  et  al.  (1973). 

Both  the  P time  term  and  the  PMP  time  showed  crustal  thickening 
under  the  LASA  center  having  a synclinal  form  and  trending  in 
the  E-NE  direction,  although  the  crustal  thickness  estimated 
from  the  Pn  data  was  not  contour able  because  of  large  fluctuations. 

The  P P results  gave  2 to  4 km  crustal  thinning  toward  north  and 

M 

4 to  6 km  thinning  toward  south  from  the  array  center  within 
the  area  under  our  reduced  LASA  aperture  (F-ring  omitted) . 

Other  models,  less  conventional  than  the  postulated  Moho- 
undulations  have  also  been  tested  using  LASA  data.  For  example. 
Mack  (1969)  explained  the  P wave-form  at  individual  stations  as 
a superposition  of  multiples  caused  by  multi-pathing  due  to 
inhomogeneities  under  the  array , although  a complete  physical 
structure  which  explains  the  whole  data  set  simultaneously 
was  not  forwarded.  On  the  other  hand,  using  both  amplitude 
and  phase-delay  fluctuation  at  various  frequencies,  Aki  (1973) 
and  Capon  (1974)  successfully  characterized  the  lithosphere 
under  the  LASA  as  a Chernov  (1960)  random  medium  with  correlation 
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distance  10  to  12  km  and  r.m.s.  velocity  fluctuations  of 
2 to  4%.  Similar  results  were  obtained  by  Berteussen  et  al.  (1975) 
who  also  found  that  a Chernov  random  medium  model  could 
explain  about  60%  of  the  variance  of  the  observed  amplitude 
fluctuations.  These  results  are  in  good  harmony  with  our 
estimate  of  the  lower  limit  of  true  velocity  fluctuation  of 
3.2%  obtained  in  the  preceding  section. 

Let  us  now  describe  our  result  using  the  dual  presentation 

of  the  generalized  inverse  and  stochastic  inverse.  Both 
solutions  for  Layer  1,  the  upper  crust,  show  a low-velocity 
anomaly  (Pigs.  4 and  9)  to  the  north  of  the  center  subarray 
AO  and  trending  roughly  N60'E.  It  coincides  well  with  the 
low  gravity  anomaly  shown  in  Fig.  2.  There  is  no  known 
geologic  feature  corresponding  to  this  anomaly.  Since  the 
sediments  are  nearly  uniform  in  this  area  (Brown  and  Poort, 

1965) , the  low  velocity  anomaly  must  be  confined  within  the 
basement  rock.  On  the  other  hand,  the  high  velocity  anomaly 
to  the  west  may  be  associated  with  a surface  feature,  i.e. 
the  Porcupine  dome. 

The  anomalies  of  Layer  2,  the  lower  crust,  are  shown  in 
Figs.  5 and  10.  A tone  of  low-velocity  anomaly  trending  N60“E 
is  again  clearly  shown  in  the  stochastic  inverse  solution. 

The  same  trend  direction  is  manifested  in  the  generalized  inverse 
solution  by  a zone  of  significantly  high  velocity  blocks.  We 
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notice  that  the  magnitude  of  velocity  perturbation  for  the 
generalized  inverse  solution  is  remarkably  small  except  for 
peripheral  blocks.  It  is  interesting  to  recall  that  the  lower 
crust  under  NORSAR  was  also  smooth  as  compared  to  the  upper 
crust  or  the  uppermost  mantle  (Paper  1) . 

The  N6Q°E  trend  of  anomaly  pattern  in  the  crust  persists 
to  the  upper  mantle  to  the  depth  of  110  km  as  shown  in  the 
velocity  anomalies  of  Layer  3 and  4 (Figs.  6,  7,  11*  and  12). 

In  all  these  solutions*  the  southeast  part  of  the  area  shows 
high  velocity.  The  boundary  between  the  zone  of  significantly 
high  and  the  ione  of  significantly  low  blocks  in  the  generalized 
inverse  solution  (Figs.  6 and  7)  agrees  well  with  the  E to  NE 
striking  zero  contours  of  the  stochastic  inverse  solution 
(Figs.  11  and  12).  Differences  between  the  two  solutions  are 
most  pronounced  in  the  NW  corner  of  the  area  sampled  where  the 
generalized  inverse  solution  continues  to  be  low  while  the 
stochastic  inverse  solution  exhibits  a moderately  high-velocity 

< 

anomaly  for  both  Layer  3 and  4.  The  anomaly  contour  pattern 
obtained  for  the  stochastic  inverse  is  very  similar  to  the 
map  of  crustal  thickness  given  by  Iyer  and  Healy  (.1972) , 

The  generalized  inverse  and  stochastic  inverse  gave 
essentially  the  same  results  for  Layer  5 except  for  the  NW 
corner.  In  addition  to  the  E to  NE  anomaly  trend  still  per- 
sisting in  the  eastern  part  of  the  LASA  siting  area*  we  find 
a secondary  trend  to  the 
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N to  NW  in  the  western  area.  This  pattern  may  be  recognized, 
though  vaguely,  in  Iyer  and  Healy's  "second  order  residual" 

projected  to  a depth  of  125  km. 

The  most  cor  spicuous  feature  of  the  estimated  velocity 
anomalies  under  LASA  is  the  N60°E  trend  which  persists  from 
the  upper  crust  to  a depth  greater  than  100  km.  This  result 
cannot  be  due  to  the  vertical  smoothing,  because  the  inter- 
layer resolution  was  perfect  for  the  generalized  inverse, 
and  very  good  for  the  stochastic  inverse.  In  order  to  obtain 
a coherent  vertical  image  of  the  anomalies  presented,  vertical 
cross-sections  of  the  two  solutions  were  constructed  along 
a profile,  indicated  by  a dashed  line  in  Fig.  2,  roughly 
perpendicular  to  the  dominant  trend  direction  N60°E  and  is 
depicted  in  Figs.  15  and  16.  In  both  cases,  we  see  high 
velocity  to  SE  and  low  velocity  to  NW  with  the  boundary 
steeply  dipping  toward  SE.  In  the  stochastic  inverse,  it 
looks  as  if  a sheet  of  low  velocity  zone  is  sandwiched  between 
two  high-velocity  regions.  On  the  other  hand,  the  low  velocity 
zone  of  the  generalized  inverse  solution  is  horizontally 
broader,  and  does  not  appear  to  be  a continuous  sheet  across 

* 

the  Moho. 

In  any  case,  it  is  well  established  that  the  same  trend 
in  the  velocity  anomalies  persists  over  nearly  the  entire 
thickness  of  the  lithosphere  with  the  lateral  scale-length  of 
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50  to  70  km.  The  maximum  velocity  contrast  is  about  4% 
for  the  stochastic  inverse  solution  and  about  10%  for  the 
generalized  inverse.  Considering  the  smoothing  effect  of  the 
former,  and  larger  random  errors  of  the  latter,  the  true 
velocity  contrast  probably  lies  between  the  above  two  values. 

What  is  the  cause  of  this  giant  anomaly  trending  N60°E? 
The  change  in  Moho  elevation  proposed  in  earlier  works  can 
at  best  account  for  the  anomalies  in  Layer  2 and  3.  The  Moho 
topography  obtained  in  earlier  works  is  likely  to  be  over- 
estimated, because  their  model  formulation  projects  the  effect 
of  inhomogeneities  other  than  Moho  topography  into  the  final 
solution.  This  will  introduce  a strong  bias  because  the 
upper-crust  and  the  lower  lithosphere  seem  to  share  the  same 
anomaly  pattern  as  the  proposed  Moho  topography.  The  sharp 
change  in  crustal  thickness  (10  km  within  the  horizontal 
distance  of  50  km)  proposed  in  earlier  studies  is  indeed 
drastic.  It  is  also  very  puzzling,  because  it  is  net 
reflected  at  all  in  the  surface  geology  nor  in  the  gravity 
anomalies.  Whatever  causes  the  anomaly,  the  trend  of  N60°E 
should  give  the  clue.  In  Fig.  2 we  notice  that  the  trends  of 
Weldon  and  Brockton  faults  are  close  to  N60°E  (Brown  and 
Poort,  1965),  and  according  to  Stone  (1965)  represents  a 
dextral  wrench  fault  belonging  to  the  system  of  wrench  faults 
prevailing  throughout  the  Rocky  Mountains  being  caused  by 
compression  in  the  NEE-SWW  direction.  Brinkworth  and 
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Kleinkopf  (1972)  show  some  correlation  between  the  strike 
of  wrench  faults  and  the  trend  of  Bouguer  gravity  anomaly 
as  shown  in  Fig.  2.  Since  the  NEE-SWW  compression  is 
consistent  with  the  two  major  orogenies  in  the  western  U.S.; 
the  Nevadan  orogeny  (Late  Jurassic-Early  Cretaceous)  and  the 
Laramide  orogeny  (Upper  Cretaceous  to  Tertiary)  we  are 
tempted  to  associate  the  N60°E  trending  velocity  anomaly 
with  shearing  caused  by  the  orogenies.  The  shear  may  be 
concentrated  in  a zone  of  weakness  which  in  turn  causes  a 
broad  low  velocity  anomaly  throughout  the  lithosphere. 


Summary 

We  have  used  the  novel  inversion  approach  of  Aki  et  al . , 
1975  in  combination  with  LASA  P wave  residuals  for  teleseismic 
events  to  estimate  the  three-dimensional  seismic  velocity 
anomalies  in  the  lithosphere  to  a depth  of  140  km  beneath 
the  array.  The  fit  to  the  observations  is  excellent  as 
indicated  by  the  final  residuals  being  comparable  to  the 
measurement  errors.  Preference  was  given  to  a dual  presenta- 
tion of  the  generalized  and  stochastic  inverse  solution  for 
the  estimated  velocity  anomalies,  because  this  facilitates 
the  geological  interpretation,  and  at  the  same  time  gives 
the  reader  greater  freedom  in  judging  the  results  than  a 
single  'optimal'  solution.  The  most  conspicuous  feature  of 
the  seismic  structure  under  LASA  is  a low  velocity  anomaly 
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in  the  central  and  NE  part  of  the  array  siting  area  with 
a N60°E  trend  and  persisting  from  the  upper  crust  to  depths 
greater  than  100  km.  We  are  tempted  to  associate  this 
low  velocity  zone  with  wrench  faulting  and  shearing  during 
the  Nevadan  and  Laramide  orogenies.  The  sharp  change  in 
crustal  thickness  proposed  in  earlier  studies  is  likely 
to  be  significantly  overestimated  because  the  upper  crust 
and  the  lover  lithosphere  appear  to  share  the  same  anomaly 
pattern  as  the  proposed  Moho  topography. 

Our  results  also  show  that  the  velocity  fluctuations  are 
of  the  same  order  in  the  crust  as  in  the  upper  mantle 

to  a depth  of  140  km.  The  root-mean- square  of  the  true 
velocity  fluctuations  in  the  lithosphere  under  LASA  is  at 
least  3.2%  and  thus  in  agreement  with  similar  estimates  of  ca 
2%  based  on  the  Chernov  random  medium  model. 
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Table  3e.  Layer  5 tabulation  for  estimated  generalised  inverse  velocity 
perturbations.  Otherwise  caption  as  for  Tables  3a  and  3b. 
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Table  5a.  Layer  1 tabulation  of  diagonal  elements  for  the  resolution 
matrix  for  the  stochastic  inverse  solution.  The  upper  left 
block  corresponds  to  the  NW  corner  of  the  grid,  while  the 
center  subarray  block  is  squared.  Initial  model  parameters 
for  Layer  1 are  listed  in  Table  2. 
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Table  5d.  Layer  4 tabulation  of  diagonal  elements  for  the  resolution 
matrix  for  the  stochastic  inverse  solution.  Caption  as  for 
Table  5a. 
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Figure  Captions 

Figure  1.  The  initial  layered  model  with  constant  thickness 
and  average  velocity  vQi  is  divided  into  many  blocks. 

Our  problem  is  to  estimate  the  velocity  perturbation  jn 
for  each  block  using  the  P travel  time  data  for  many 
teleseismic  events  observed  at  the  receivers  above  the 
model.  The  velocity  perturbation  _m  for  a particular 
block  in  a particular  layer  represents  an  average  over 
all  the  rays  having  most  of  the  unperturbed  path  in 
the  layer  within  the  block. 

Figure  2.  Map  of  LASA  subarray  centers  and  Bouguer  gravity 

anomalies.  The  dashed  line  marks  the  profile  along  which 
vertical  cross-sections  of  velocity  anomalies  are  shown 
in  Figures  15  and  16. 

Figure  3.  Epicenter  locations  of  the  178  events  used  in  analysis. 

Figure  4.  Generalized  inverse  solution  for  Layer  1.  The 

numbers  show  the  fractional  velocity  perturbation  in  per 
cent  of  the  average  layer  velocity  given  in  Table  2.  The 

n 

snaded  areas  correspond  to  the  magnitude  of  solution  greater 
than  twice  the  standard  error  which  is  listed  in  Table  3. 

The  letters  L and  H refer  to  low  and  high  velocity  anomalies 
respectively.  The  central  subarray  is  marked  by  a larger 
solid  circle. 

Figure  5.  Generalized  inverse  solution  for  Layer  2.  The 

unresolved  block  is  encircled.  See  Figure  4 for  explanation 
of  symbols. 
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Figure  6.  Generalized  inverse  solution  for  Layer  3.  See 
Figures  4 and  5 for  explanation  of  symbols. 

Figure  7.  Generalized  inverse  solution  for  Layer  4,  See 
Figures  4 and  5 for  explanation  of  symbols. 

Figure  8.  Generalized  inverse  solution  for  Layer  5.  See 
Figures  4 and  5 for  explanation  of  symbols. 

Figure  9.  Stochastic  inverse  solution  for  Layer  1.  The 
numbers  show  the  fractional  velocity  perturbation  in 
per  cent  of  the  average  velocity  given  in  Table  2.  The 
values  represent  the  4 point  average  over  two  block 
configurations,  one  centered  at  the  central  subarray 
(marked  by  a large  solid  circle)  and  the  other  shifted 
by  a half-block  length  towards  southwest.  The  letters 
L and  H mark  areas  of  low  and  high  velocity  anomalies 
respectively.  The  contours  are  drawn  at  1%  intervals. 

Figure  10.  Stochastic  inverse  solution  for  Layer  2.  See 
Figure  9 for  explanation  of  symbols. 

Figure  11.  Stochastic  inverse  solution  for  Layer  3.  See 
Figure  9 for  explanation  of  symbols. 

Figure  12.  Stochastic  inverse  solution  for  Layer  4.  See 
Figure  9 for  explanation  of  symbols. 

Figure  13.  Stochastic  inverse  solution  for  Layer  5.  See 
Figure  9 for  explanation  of  symbols. 
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Figure  14.  Selected  elements  of  the  resolution  matrix  in  % 
for  the  stochastic  inverse  solution.  The  diagonal 
elements  are  enclosed  by  thicker  lines.  For  example, 
the  five  tic-tac-toes  under  A are  the  elements  of  row 
vector  corresponding  to  the  diagonal  element  in  Layer  3 
at  geographic  location  "A"  marked  on  the  right.  See 
text  for  detailed  explanation. 

Figure  15.  Generalized  inverse  solution  for  the  vertical 

cross-Section  along  the  profile  indicated  in  Figure  2. 
The  numbers  were  obtained  by  interpolation  from  those 
given  in  Figures  4 through  8.  The  letters  H and  L refer 
to  high  and  low  velocity  anomalies  respectively. 

Figure  16.  Stochastic  inverse  solution  for  the  vertical 

cross-section  along  the  profile  indicated  in  Figure  2. 
The  numbers  are  directly  obtained  from  Figures  9 through 
14.  The  letters  H and  L refer  to  high  and  low  velocity 
anomalies  respectively. 


-150- 


t 


LAYER  3 (50-80  KM) 


LASA 


FIGr  II 
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PARTS  OF  RESOLUTION  MATRIX  NEAR  DIAGONAL  ELEMENTS 
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4 . 4 Elevation  of  the  Olivine-Spinel  Phase  Change  in  Su bduc t e d 

Lithosphere:  Seismic  Evidence  by  S.C.  Solomon  and 

K.T.  Paw  U 
ABSTRACT 

The  top  of  the  olivine-spinel  phase  change  in  subducted 
oceanic  lithosphere  can  be  located  by  the  travel  times  of 
seismic  waves  which  have  propagated  through  the  slab.  P-wave 
travel  time  residuals  from  deep  earthquakes  in  the  Tonga 
island  arc  observed  at  Australian  seismic  stations  are  grouped 
according  to  the  depth  of  the  earthquake.  The  change  in  mean 
residuals  with  a change  in  earthquake  depth  is  related  to  the 
velocity  contrast  between  slab  and  normal  mantle  at  that 
depth.  The  curve  mean  residual  versus  earthquake  depth 
displays  a region  of  markedly  increased  slope  between  earth- 
quake depths  of  about  250  and  350  km.  The  most  probable 
explanation  of  this  observation  is  an  elevation  by  100  km  of 
the  olivine-spinel  phase  change  within  the  relatively  cooler 
slab.  No  evidence  was  found  for  vertical  displacements 
within  the  slab  of  any  deeper  phase  changes. 

A temperature  contrast  between  slab  and  normal  mantle 
of  about  1000°C  at  250  km  depth  is  implied.  This  finding 
confirms  current  thermal  models  for  subducted  lithosphere 
but  is  inconsistent  with  the  global  intraplate  stress  field 
unless  only  a few  percent  of  the  negative  buoyancy  force  at 
subduction  zones  is  transmitted  to  the  surface  plates. 


— 3 6 f> -• 


4.5  P Wave  Travel  Times  from  Deep  Earthquakes  by  M . K . 


Sengupta  and  B.R.  Julian 


ABSTRACT 


A study  of  teleseismic  P wave  travel  times  has  been 


made  using  data  from  deep  (450-650  km)  earthquakes  exclusively, 


for  which  the  effect  of  near  source  upper  mantle  heterogeneity 


is  expected  to  be  small  and  location  parameters  are  considered 


reliable.  About  3,300  arrival  times  selected  from  bulletins 


and  winnowed  to  remove  discrepant  data,  were  used.  The  shape 


of  the  travel  time  curve  and  the  average  station  anomalies 


found  are  in  general  agreement  with  the  results  of  other  works 


(e.g.  Herrin  et  al . , 1968)  but  the  scatter  of  the  data  is 


only  about  half  as  large  <M>.6  tec),  verifying  that  deep 


earthquakes  are  indeed  more  suitable  than  shallow  ones  for 


travel  time  studies.  The  shape  of  the  travel  time  curve  is 


determined  within  0.1-0. 2 seconds  out  to  epicentral  distances 


of  70°,  but  the  absolute  level  of  th2  curve  remains  uncertain 


within  about  ±1  second.  This  'baseline  correction'  has  been 


estimated  from  Nevada  Test  Site  explosion  data,  for  which 


the  upper  mantle  velocity  structure  near  the  source  is  known 
and  could  be  corrected  for.  The  results  of  this  procedure 


indicate  that  the  Jef f reyr.-Bullen  (1940)  times  are  about 


1.5  seconds  late  but  the  model  of  Herrin  et  al.  (1968)  has 


approximately  the  correct  baseline  of  P travel  times  at 


least  for  focal  depth  of  550  km. 
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The  station  anomalies  of  this  study,  represented  by  a 
constant,  correlate  very  well  with  large-scale  tectonic 
features  and  also  give  a better  fit  to  the  travel  time  data 
than  has  been  achieved  in  other  studies  which  included 
azimuthal  station  terms. 

Systematic  regional  variations  in  travel  time  of  1 
second  or  more  occur  beyond  70°,  and  indicate  that  signifi- 
cant lateral  variation  occurs  in  the  lower  mantle  (depth  > 

2000  km) . A conservative  lower  bound  of  about  1%  can  be 
placed  on  the  magnitude  of  these  variations  in  compressional 
velocity. 

4 . 6 Three  Dimensional  Model  of  Seismic  Velocity  Variation 
in  the  Earth's  Mantle  by  M.K.  Sengupta  and  M.N.  Toksflz 
Three  dimensional  velocity  models  for  the  earth's  mantle 
were  obtained  satisfying  P,  S and  some  PcP  and  ScS  travel 
time  anomalies  from  deep  focus  earthquakes.  A method  of 
successive  approximation  was  used  to  compute  the  uniform 
relative  velocity  perturbation  over  three  dimensional  blocks 
of  size  10°xl0°  in  longitude  and  latitude,  and  500  km  in 
thickness.  Features  of  these  velocity  perturbations  indicated 
that  lateral  heterogeneity  is  the  most  pronounced  in  the  upper 
mantle  and  near  the  core-mantle  boundary.  The  upper  mantle 
anomalies  are  correlated  with  surface  tectonic  features  but 
these  are  not  correlated  with  deep  mantle  anomalies. 
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Global  and  regional  analyses  of  geophysical  data  indicate 
that  there  are  lateral  heterogeneities  in  the  earth's  upper 
mantle  as  well  as  in  the  lower  mantle  (Toksfiz  et  al^  , 1967; 
Davies  and  Sheppard,  1972;  Julian  and  Sengupta,  1973;  Jordan 
and  Lynn,  1974;  Wright  and  Lyons,  1975).  In  this  paper  we 
present  three  dimensional  seismic  velocity  models  of  the 
earth's  mantle  derived  from  good  quality  P,  S and  some  PcP 
and  ScS  travel  time  data  from  deep  focus  earthquakes.  In 
earlier  papers  Davies  and  Sheppard  (1972)  and  Julian  and  Sen 
Gupta  (1973) , using  array  diagram  and  travel  time  anomalies 
along  different  ray  bundles,  concluded  that  the  lowermost 
mantle  (.2000-3000  km)  is  comparatively  more  heterogeneous 
laterally  than  the  middle  mantle  (1000-2000  km) . In  the  present 
paper,  we  look  at  the  inverted  model  for  the  whole  mantle 
which  was  divided  into  blocks.  Since  velocity  variations  are 
generally  small  we  determine  percentage  perturbation  of 
velocities  relative  to  an  average  model. 


Data 


Travel  times  of  compressional  (P) , shear  (S) , and  some 
core-reflected  phases  (PcP,  ScS)  data  were  read  from  World 
Wide  Seismic  Station  Network  (WWSSN)  seismograms  of  12  deep 
focus  (depth  > 500  km)  earthquakes  under  the  Indonesian  arc, 

Solomon  Islands,  New  Hebrides,  Kuril  arc,  Tonga  arc,  and 

* 

western  coast  of  South  America.  To  expand  the  global  coverage, 
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we  also  included  data  from  two  intermediate  depth  earthquakes 
from  Greece  and  Hindikush.  To  increase  the  ray  path  distri- 
butions, we  selected  P wave  arrivel  time  data  from  international 
bulletins  for  ether  stations  around  the  world.  All  of  the 

r 

travel  time  data  from  these  14  earthquakes  were  assumed  to 
be  free  from  any  severe  source  bias  as  these  earthquakes  were 
among  the  deepest  ones  in  the  respective  Benioff  zones. 

The  earthquakes  were  relocated  after  constraining  the  depths 
from  readings  of  the  depth  phases  (pP,  sP)  . 


Method  of  Inversion 

The  observed  travel  times,  averaged  over  2°  distance 
cells,  were  fitted  first  with  an  average  radial  velocity  model. 
Then,  the  difference  between  the  observed  and  calculated 
travel  times  was  averaged  for  each  station  to  obtain  the 
"station  anomalies".  The  "station  anomalies"  incorporate 
the  crustal  anomalies  and  the  short  wavelength  lateral  vari- 
ations in  the  upper  mantle.  The  remaining  travel  time  anomalies 
were  fitted  through  constant  velocity  perturbations  in  three 
dimensional  blocks  of  sizes  10°  in  longitude,  10°  in  latitude, 
and  500  km  in  thickness.  About  2000  such  blocks  were  sampled 
by  1490  P and  PcP  travel  time  data  and  1400  blocks  by  314  S 
and  ScS  travel  time  data. 

To  compute  the  velocity  perturbation  in  each  block,  we 
used  a method  of  successive  approximation.  In  this  method, 
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the  rth  approximation  of  the  velocity  perturbation  for  the 
mth  block  (=  pm)  is  given  by  the  weighted  least  square  solu- 
tion of  the  system  of  n linear  equations  of  the  form: 


AT??S  + cf,  = KL_  P, 
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ID 
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ijm  ^m 


(1) 


where  n is  the  number  of  ray  paths  passing  through  the  mth 
block,  AT?jS  is  the  observed  travel  time  anomaly  for  ith 
earthquake  and  jth  station,  c^j  is  the  rth  approximation  of 
the  correction  to  AT?^S  for  the  perturbation  of  velocity  in 
other  blocks,  K^jm  is  the  rth  approximation  of  "data  kernel" 
of  the  mth  block  corresponding  to  AT?jS.  These  data  kernels 
are  actually  the  calculated  travel  times  of  the  corresponding 
rays  in  the  given  block.  During  the  computation  of  velocity 
perturbation,  weights  were  given  to  each  datum  according  to 
the  fraction  of  the  total  travel  time  the  ray  spends  in  the 
given  block.  The  computed  relative  velocity  perturbation  was 
also  constrained  not  to  exceed  5%  (or  1%  if  there  is  only  one 
datum)  in  any  individual  cycle  of  iteration.  Altogether 
three  cycles  were  performed  before  convergence  was  reached. 

The  model  in  our  method  of  inversion  depends  on  the 

sequence  in  which  the  respective  blocks  were  chosen  during 

inversion.  So,  three  different  sequences  were  chosen.  In 

the  first  trial  (MODEL  1) , the  blocks  were  sequenced  in  an 

ordered  way  starting  from  the  core-mantle  boundary  and  ter- 

* 

minating  on  the  free  surface;  in  the  second  trial  (MODEL  2) 


n 

-165- 

the  above  sequence  was  reversed.  Finally  in  MODEL  3,  the 

blocks  were  sequenced  randomly . After  fit  through  these  three 

dimensional  models,  r.m.s.  deviation  of  1490  P-wave  travel 

time  residuals  was  about  0.4  seconds,  comparable  to  the 

€ 

standard  deviation  of  the  reading  error.  Since  S-wave  travel 
time  data  are  fewer  than  those  of  P-waves  (314  vs.  1490), 
the  S-wave  data  were  over fitted  by  the  models  and  produced 
r.m.s.  deviation  of  only  0.2  seconds. 

Result  of  Inversion 

As  the  absolute  magnitude  of  the  velocity  perturbation 
depends  on  the  size  and. sampling  of  the  blocks,  we  discuss 
only  the  "relative  variation''  of  the  velocity  perturbations 
with  depth  (Fig.  1)  and  with  latitude  and  longitude  (Fig.  2). 

I 

The  velocity  perturbations  for  the  surface  layer  (0-500  km) 
shown  in  these  two  figures  include  a contribution  corresponding 
to  the  average  of  station  anomalies  in  each  block.  All  three 
models  shown  in  Fig.  1 show  a relative  maximum  for  the 
surface  layer  (0-500  km)  and  for  the  layer  near  the  core- 
mantle boundary  indicating  that  these  two  layers  perhaps  have 
the  largest  lateral  heterogeneity  of  compressional  waves. 

For  shear  waves,  the  greatest  heterogeneity  in  the  upper 
mantle  (0-500  km)  is  verified  by  all  models.  Heterogeneity 
near  the  core-mantle  boundary  is  not  as  clearly  defined. 

The  compressional  wave  data  have  sufficiently  broad 
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distribution  to  investigate  the  lateral  variations  of  the 
velocities.  In  Fig. 2 we  show  the  relative  P-wave  velocity 
perturbations  (in  percent)  in  the  upper  mantle  and  near  the 
core-mantle  boundary.  These  velocity  perturbations  were 
obtained  using  inversion  with  random  sequencing  of  blocks 
(MODEL  3) . This  figure  illustrates  that  the  pattern  of  P 
velocity  anomalies  for  the  upper  mantle  of  our  model  is  in 
agreement  with  other  seismological  data  and  tectonic  features. 

Generally  the  tectonically  active  regions  are  character- 
ized by  slower  average  velocities  (negative  perturbations) 
and  shield  arec^s  by  higher  velocities  (positive  perturbations) 
in  the  outer  500  km.  This  is  obvious  in  North  America  with 
a contrast  between  the  eastern  and  western  halves.  In  Africa 
and  in  Europe,  shield  areas  stand  out  as  regions  of  higher 
average  velocities.  The  contrasts  between  the  stable  Indian 
shield  and  the  Himalayas  and  Tibet  and  Southeast  Asia  are 
similar  to  those  observed  from  surface  wave  studies  (Bird  and 
Toksfiz,  1975) . In  Japan  our  results  are  similar  to  those 
obtained  by  Kanamori  (1970)  from  regional  studies  and  in 
Australia  to  those  obtained  by  Cleary  (1967).  Since  the  top 
(0-500  km)  layer  is  sampled  directly  under  the  stations,  the 
data  are  co.ifined  to  continental  areas  and  to  a few  islands 
where  seismic  stations  exist.  Under  Hawaii  the  average 
velocities  are  low.  In  oceanic  areas  there  are  no  data. 

For  the  core-mantle  boundary  region,  it  is  difficult  to 


pH 


;• 


verify  any  model  because  of  both  scarcity  of  proper  geophysical 

data  and  of  ambiguity  in  interpretation  of  those  data.  In 

1 

our  P velocity  model  near  the  core-mantle  boundary,  the  central 

Pacific  region  stands  out  as  having  small  variation  in  velocity 

• *m 

(within  ±5%)  compared  to  our  spherically  symmetric  model. 


Discussion  and  Conclusions 

Inverted  velocity  models  are  dependent  on  the  size  of 
the  blocks  and  the  sampling,  or  the  number  of  rays  passing 
through  each  block.  Block  sizes  must  be  chosen  small  enough 
to  delineate  anomalies,  and  large  enough  to  have  a sufficient 
number  of  rays  passing  through  them.  The  bl^ck  sizes  chosen 
in  this  study  (of  the  order  of  500  km)  are  generally  larger 
than  scale  lengths  of  mantle  heterogeneities.  Detailed 
studies  of  subduction  zones  (Toksfiz  et  aa. , 1971)  and  the 
core-mantle  boundary  (Doornbos,  1975)  indicate  that  scale 
length  of  heterogeneities  in  the  earth's  mantle  is  between 
10-50  km.  Thus  our  models  grossly  average  over  such  smaller 
heterogeneities.  Actual  local  perturbations  could  be  much 
larger  than  those  given  in  Figs.  1 and  2. 

During  the  inversion,  we  found  that  the  number  of  rays 
sampling  each  block  also  had  an  influence  on  the  results. 
Blocks  with  few  hits  (<5)  generally  had  two-fold  larger 
velocity  perturbations  than  the  blocks  sampled  by  more  data. 
This  is  partly  due  to  the  artifact  of  inversion  especially 
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for  the  blocks  with  only  one  hit.  Another  reason  for  it  nay 
lie  in  the  largeness  of  the  block  size  compared  to  the  scale 
length  of  heterogeneity  in  the  real  earth  as  mentioned  in 
the  previous  paragraph.  Deleting  the  blocks  with  few  hits 

r'* 

(<5)  did  not  alter  the  general  picture.  The  relative  varia- 
tion of  P wave  velocity  perturbations  with  depth  followed  a 
similar  trend  as  shown  in  Fig.  1 (top) . 

The  basic  conclusions  that  are  consistent  with  the 
analysis  of  our  travel  time  data  are: 

(1)  Lateral  heterogeneity  is  most  pronounced  in  the 
upper  mantle  and  near  the  core-mantle  boundary  of  the  earth; 

(2)  In  the  upper  mantle  the  velocity  variations  are 
correlated  with  tectonic  features; 

(3)  Comparing  the  features  of  the  P velocity  model  for 
the  upper  mantle  with  those  near  the  core-mantle  boundary, 
one  finds  no  relationship  between  the  surface  tectonics  and 
the  deep  mantle  structure. 
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Figure  Captions 

Figure  1.  R.m.s.  relative  perturbations  Cin  percent)  of  P 
and  S wave  velocity  for  each  500  km  depth  interval. 

I I 

Figure  2.  Three  dimensional  P wave  velocity  models  for  the 
. upper  mantle  and  for  near  the  core-mantle  boundary. 

Each  color  represents  a range  in  relative  velocity 
perturbation  (in  percent) . The  positive  values  in 
relative  velocity  perturbation  correspond  to  higher 
velocity  than  spherically  symmetric  model  and  the 
negative  perturbation  to  lower  velocity. 
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SCATTERING  AND  ATTENUATION 


5.1  Origin  of  Coda  Waves:  Source,  Attenuation,  and 

Scattering  Effects  by  K.  Aki  and  B.  Chouet 
Coda  waves  from  small  local  earthquakes  are  interpreted 
as  backscattering  waves  from  numerous  heterogeneities  dis- 
tributed uniformly  in  the  earth's  crust.  Two  extreme  models 
of  the  wave  medium  that  account  for  the  observations  on  the 
coda  are  proposed.  In  the  single  backscattering  model  the 
scattering  is  considered  to  be  a weak  process,  and  the  loss 
of  seismic  energy  by  scattering  is  neglected.  In  the  dif- 
fusion model  the  seismic  energy  transfer  is  considered  as  a 
diffusion  process.  Both  models  lead  to  similar  formulas 
that  allow  an  accurate  separation  of  the  effect  of  earthquake 
source  from  the  effects  of  scattering  and  attenuation  on  the 
coda  spectra.  A unique  difference  was  found  in  the  scaling 
law  of  earthquake  source  spectra  between  central  California 
and  western  Japan,  which  may  be  attributed  to  the  difference 
in  inhomogeneity  length  of  the  earth's  crust.  The  Q of 
coda  waves  in  the  two  regions  is  strongly  frequency  dependent 
with  values  increasing  from  50-200  at  1 Hz  to  about  1000- 
2000  at  20  Hz.  This  observation  is  interpreted  as  a combined 
effect  of  variation  of  Q with  depth  and  frequency-dependent 
composition  of  coda  waves  described  below.  The  turbidity 
coefficient  of  the  lithosphere  required  at  1 Hz  to  explain 
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the  observed  coda  as  body  wave  scattering  is  orders  of  magni- 
tude greater  than  previously  known  values  such  as  those 
obtained  by  Aki  (1973)  and  Capon  (1974)  under  the  Montana 
LASA  from  the  amplitude  and  phase  fluctuations  of  teleseismic 
P waves.  From  the  high  attenuation  and  turbidity  obtained 
at  this  frequency  we  conclude  that  at  around  1 Hz  tfr e coda 
is  made  of  backscattering  surface  waves  from  heterogeneities 
in  the  shallow,  low-Q  lithosphere.  The  high  Q observed  for 
the  coda  at  frequencies  higher  than  10  Hz,  on  the  other 
hand  eliminates  the  possibility  that  these  waves  are  back- 
soattering  surface  waves.  We  conclude  that  at  these  high 
frequencies  the  coda  must  be  made  of  backscattering  body 
waves  from  heterogeneities  in  the  deep  lithosphere.  The 
low  turbidities  found  for  deep  earthquake  sources  under 
western  Japan  are  consistent  with  this  model  of  coda  wave 
generation. 


5 . 2 Synthesis  of  Teleseismic  P-waves  from  Sources  near 

Sinking  Lithospheric  Slabs  by  R.W,  Ward  and  K.  Aki 

■ 

ABSTRACT 

A wave  theory  method  is  used  to  determine  the  effect,  of 
a sinking  lithospheric  slab  on  short-period  and  long-period 
waves.  We  consider  a simplified  model  of  the  lithospheric 
flab  with  ? 10%  velocity  contrast  and  compute  both  short- 
period  and  long-period  theoretical  seismograms  from  a P-wave 
source  located  in  or  near  the  slab.  For  this  model,  the 
ray-theoretical  amplitude  agrees  quite  well  wit  the  short- 
period  amplitude.  In  the  ray-theoretical  shadow  zone  the 
long-period  seismograms  typically  have  amplitudes  50%  (or 
greater)  of  the  direct  P-wave  fmplitude  and  exhibit  waveform 
broadening.  Similar  waveform  broadening  has  been  attributed 
to  the  dynamics  of  earthquake  faulting.  The  effect  of  the 
lithosphere  on  long-period  waves  from  nearby  sources  must  be 
taken  into  account  in  studies  which  utilize  the  observed 
variation  in  waveform  broadening  to  infer  earthquake  source 
dynamics . 
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5 . 3 Inversion  Schemes  for  Surface  Wave  Attenuation  and  Q 

in  the  Crust  and  the  Mantle  by  W.B.  Lee  and  S.C.  Solomon 
Seismic  wave  attenuation,  or  t)  \ is  one  of  the  para- 
meters traditionally  used  to  distinguish  ..ithosphere  from 
asthenosphere.  As  a means  of  inverting  surface  wave  attenu- 
ation measurements  to  obtain  Q ^ as  a function  of  depth,  a 
classical  linear  inverse  problem,  we  introduce  a set  theoretical 
approach.  The  set  theoretical  approach  includes  both  the 
square  matrix  inverse  and  the  linear  programming  method,  which 
agree  in  concept  but  differ  in  that  the  constraints  consist 
of  mean  or  extreme  hyperplanes,  respectively.  The  set 
theoretical  approach  to  inverting  attenuation  observations 
has  two  advantages:  (1)  The  lithosphere  thickness  may  be 

readily  estimated  using  a geometric  visualization  of  the 
square  matrix  inverse,  which  identifies  the  feasibility  of 
a solution.  (2)  The  approach  allows  more  strict  regulation 
of  solution  parameters  than  does  the  least  square  inverse, 
which  often  leads  to  a negative  solution  for  inaccurate  and 
sparse  attenuation  data.  We  apply  the  set  theoretical 
approach  to  inverting  the  measured  attenuation  of  Love  and 
Rayleigh  waves  in  western  and  east-central  United  States. 
Identifying  the  base  of  the  lithosphere  as  the  depth  at 
which  Q increases  sharply,  we  obtain  lithosphere  thick- 
nesses of  80+20  km  and  130+30  km  in  the  two  regions, 
respectively.  Q-^  in  the  asthenosphere  is  about  a factor 
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